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A FINITE DIFFERENCE METHOD FOR THE SOLUTION 
OF THE TRANSONIC FLOW AROUND HARMONICALLY 
OSCILLATING WINGS 

by F. Edward Ehlers 
Boeing Commercial Airplane Company 


SUMMARY 


A finite difference method is presented for the solution of the unsteady pressure distribution 
on harmonically oscillating wings in transonic flow. The differential equation for the velocity 
potential of compressible flow was simplified first by the assumption of small perturbation from a 
nearly sonic parallel flow and then by representing the flow as a sum of two separate potentials for 
the steady and unsteady flow. 

On the assumption that the amplitude of the unsteady oscillation is small, the unsteady 
transonic flow differential equation was linearized. The resulting differential equation contains 
coefficients which depend upon the mean steady flow and which must be computed beforehand 
from the solution of the nonlinear transonic small perturbation differential equation. When time is 
eliminated from the linear unsteady flow equation by the assumption of harmonic motion, then the 
resulting linear differential equation is of mixed type, being elliptic or hyperbolic in regions of the 
flow where the steady flow equation is elliptic or hyperbolic, respectively. 

Because of the similarity of the unsteady linear differential equation to the steady nonlinear 
transonic small perturbation differential equation, a technique of column relaxation similar to that 
of Murman and Cole (ref. 2) and Krupp (refs. 3 and 4) was applied to find solutions to the 
differential equation expressed in difference form. For subsonic mesh points in the flow, central 
differences were employed for all derivatives; for supersonic points, backward (or upstream) 
differences were used for derivatives in the streamwise direction only. The boundary conditions for 
the wing were applied in the plane of the projection of the wing surface and incorporated in the 
difference equations for mesh points adjacent to the oscillating wing. On the plane vortex sheet 
downstream of the trailing edge, the condition of continuity of pressure and the Kutta condition on 
the trailing edge specify relations between the jump in potential along lines in the free-stream 
direction and the value of this jump just downstream of the trailing edge. 

Integral relations based on the wing and wake boundary conditions were derived for the wing 
and the airfoil oscillating in harmonic motion in the manner of Klunker (ref. 8). These relations are 



used for the mesh far-field boundary conditions, which are updated during the computations. The 
use of far-field boundary conditions makes it possible to reduce the region of calculations, thereby 
saving computer time and storage. 

Expanding the terms in the difference equation in Taylor series about the central point and 
retaining the lowest order terms in the mesh size lead to the original differential equation plus 
additional higher derivative terms which vanish as the mesh size goes to zero. For equally spaced 
points and subsonic flow, the resulting fourth-order differential operator is of second order in the 
mesh spacing and elliptic in form. Hence, there is no dependence upon the ratio of the mesh sizes 
for stability. For locally supersonic points, the additional truncation terms are of first order in the 
mesh spacing. This term is a third derivative in the free-stream direction and introduces a simulated 
viscosity proportional to the mesh size which was found to stabilize the calculations for shocks 
developing in the flow in solving for the steady flow equations by Murman and Cole (ref. 2) and 
Krupp (refs. 3 and 4). From this study, the differencing method would appear to be stable. 
Subsequent experience has shown this to be true for the oscillating airfoil. 

The mathematical development of the method described in the foregoing for oscillating 
rectangular wings, swept wings, and airfoils is presented here in sufficient detail to be coded for 
actual numerical calculation. To test the feasibility of the method, the equations for the 
harmonically oscillating flap on a NACA 64A006 airfoil were coded for the CDC 6600. For the 
basic steady flow, the program (TEA-330) developed by Krupp in reference 4 for the transonic flow 
over a lifting airfoil was used. 

Solutions were obtained for the unsteady flow over a flat plate with an oscillating flap of a 
quarter-chord length at a free-stream Mach number of 0.8 and reduced frequencies of to = 0.06 and 
0.1794. These values correspond to frequencies of 30 and 90 cycles per second for the wind tunnel 
and model conditions in the experiments of Tijdeman and Schippers (ref. 9). The real part of the jump 
in pressure coefficient from the difference method agrees very closely with solutions from a 
subsonic kernel function method which solves the same differential equation. The imaginary part of 
the pressure coefficient has somewhat smaller negative values near the leading edge for cu=0.06, 
although* the pattern is essentially the same as the results from the kernel function method. For 
oj = 0.1794, the imaginary parts from the two methods agree more closely than for oo= 0.06. The 
convergence was slow and required about 3000 iterations to converge to a maximum difference 
between iterations of 2 x 10"^ starting from a value of zero for the perturbation potential at all 
mesh points. See note added in proof on page 34. 

The flat plate solution for a reduced frequency of 0.06 and Mach number of 0.8 was used as a 
starting solution for the unsteady flow around a NACA 64A006 airfoil with a quarter-chord 
oscillating flap at Mach 0.794 and reduced frequency c o = 0.064. Similar convergence of 2 x 10"** 
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was attained in 1200 iterations. The pattern of the pressure jump across the airfoil agrees with the 
experiments of Tijdeman and Schippers (ref. 9), although the amplitude is considerably higher. The 
solution for M = 0.804 and u) = 0.253 (corresponding to a frequency of 120 cps) about the NACA 
64A006 airfoil was obtained from the solution for M = 0.794 and a? =0.064 as starting values. 
Solutions were also found for M = 0.858 and w = 0.179 and for M = 0.853 and co = 0.060. The 
patterns for the jump in pressure coefficient from the calculations were similar to the experimental 
results of Tijdeman and Schippers in reference 9 but the amplitudes were generally higher. 

Although only a few computed results were 'obtained, the feasibility of this approach for 
finding the unsteady pressure distributions around an oscillating wing in the transonic range is 
demonstrated. This method predicts the observed pressure patterns, while linearized subsonic 
theory does not. The numerical procedure for the unsteady flow is stable and convergent even when 
the steady flow contains local supersonic regions. 

All the calculations were made with a fine fixed grid. Garabedian and Korn (ref. 1 ) found that 
first obtaining a solution of the exact steady flow potential with a coarse grid and then refining the 
grid for greater accuracy significantly reduced the number of iterations required for convergence of 
the final solution. A similar procedure should reduce the number of iterations for the unsteady 
solution as well. 

The complete matrix for the coefficients of the unknown values of the unsteady perturbation 
potential at the points of the mesh contains only five or six nonzero diagonals. The time required 
for solution also may be reduced by using a special method for solving linear equations with sparse 
matrices to solve the complete set of difference equations rather than using column relaxation. This 
may be especially useful for the coarse grid solution. 

The basic technique followed here of representing the flow as consisting of separate steady and 
unsteady flow potentials and then linearizing the unsteady differential equation by assuming small 
amplitude of oscillation can be applied directly to the complete nonlinear differential equation for 
compressible flow as well as to the nonlinear transonic small perturbation equation. The problem is 
easier for the airfoil than for the three-dimensional wing since mapping techniques such as that of 
Garabedian and Korn in reference 1 may be used to represent the boundary in a way convenient for 
the mesh system. The same finite difference scheme could be used for the unsteady flow equation 
after it is expressed in the same independent variables. The elimination of the small perturbation 
restriction for the three-dimensional wing solution is considerably more difficult, since convenient 
body coordinates are possible only for restricted classes of configurations. A suitable variable grid 
around the actual wing boundaries must be devised with appropriately defined difference equations. 
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Further computations should be made with the airfoil program in its present form. The 
sensitivity of the program to mesh size near singularities and to various ways of conducting the 
column relaxation should be studied. Further calculations are needed to gain experience in 
determining optimum values of overrelaxation for subsonic regions and underrelaxation for 
supersonic regions, and to Find criteria for the degree of convergence required for a good solution. 
Some of the discrepancy between calculated and experimental results may be due 
to the influence of the boundary layer and/or separation effects over the aft 
portion of the airfoil. Unknown problems associated with the theory or with 
the experiment, such as wind tunnel wall interference, may also be the cause 
of this discrepancy. 


INTRODUCTION 


In recent years, significant advances have been made in the theoretical treatment of 
two-dimensional steady transonic flow by such workers as Garabedian and Korn (ref. 1), Murman 
and Cole (ref. 2), Krupp (refs. 3 and 4), and Steger and Lomax (ref. 5). More recently, some 
progress has been achieved by Newman and Klunker (ref. 6), and Ballhaus and Bailey (ref. 7) 
toward the development of a useful method for predicting steady transonic airloads on finite wings. 
There are, however, no satisfactory methods for predicting transonic loads on oscillating surfaces 
for use in flutter and gust analysis of aircraft. Any completely adequate analysis for the transonic 
region should include the effects of airfoil bluntness, thickness, camber, angle of attack, and 
wing-body interference. It should include the effects of mixed subsonic and supersonic flows 
containing embedded shocks and of the boundary layer with shock wave interaction. 

A method for treatment of oscillating wings by finite difference methods is presented here. 
The assumption of small perturbations from a uniform stream near the speed of sound allows the 
reduction of the equations of motion to a simple, more tractable form and retains the necessary 
nonlinearity for describing flows with local supersonic regions. The introduction of the perturbation 
velocity potential restricts the solution to weak shocks which, for thin wings of reasonably good 
design, is not too limiting an assumption. When the flow is steady the resulting nonlinear 
differential equation reduces to the well-known transonic small perturbation differential equation 
recently studied by Murman and Cole (ref. 2) and Krupp (refs. 3 and 4). The unsteady flow 
differential equation is simplified by considering the flow as consisting of the sum of two separate 
potentials representing the steady and unsteady effects and by linearizing the differential equation 
for small amplitudes of harmonic oscillation. 
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SYMBOLS AND ABBREVIATIONS 


a,b 

c,d 

c l,dj, c 2 

c 

c 3i 

c 4i ,d 4i 

c sl ,c s2’ d sl ,d s2 

c kl> c k2’ c k3’ c k4 
D n = i3/9x n 

E 

f(x,y,t) - 
f(x,t) 

f 0 

f l 

^ij “ ^ 1 x + * ^ 1 
h 

h i 

ij,k 

*a 

^max^max’^max 

h 

Jm 


coefficients for y,z or differences corresponding to second derivatives, 
with appropriate subscripts (equation (23) in appendix) 

coefficients for x or £ difference corresponding to second derivative 
(equation (19) in appendix) 

coefficients for second-order accurate difference corresponding to first 
derivative (equations (20) and (26) in appendix) 

velocity of sound in units of Uq 

equation (92) in appendix 

equation (94) in appendix 

equation (105) in appendix 

equations (117) and ( 1 2 1 ) in appendix 

(x 1 ,x 2 ,x 3 ) = ($,T?,f) 

coefficients in difference equations, with appropriate subscripts 
instantaneous wing shape defined by Zq = 5f(x,y,t) 
instantaneous airfoil shape defined by yg = 5f(x,t) 
undisturbed wing or airfoil shape 
unsteady contribution to wing or airfoil shape 
unsteady wing boundary conditions 


z k +1 " z k 
K nv 1 K m 

x,y,z or.£,77,f subscripts for points in the mesh 
x index for first mesh point behind hinge 
maximum number of x,y,z mesh planes, respectively 
x index for trailing edge 
y mesh line just below airfoil 
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K 

k-m 

L,U 

M 

P(D) 

P(x,y) 

q 

r 

s 


S 1 


s 2 

M 0 

u 0 

u,v,w 


u 

V 

x,y,z 

x O’yo> z o 

xj.yj.zj 

a 

0 

01 

7 


(1-M 2 /M 2 e 

z mesh plane below wing 

superscripts denoting lower or upper boundary on Fjj 
free-stream Mach number 
polynomial operator, D = (Dj,D 2 ,D 3 ) 

Y>ix + icoi^i acceleration (or pressure) potential 
'o 2 /e -My- l)^0xx 



semi-chord of wing 


(yj +2 -y; +i)/h 

J m Jm 1 



time in units of s/uq 
free-stream velocity 

x,y,z velocity components of flow at beginning of the derivation of the 
equations 

K - (y + l)>Pox throughout remainder of text 
K - (y + 1 Vox + ^ throughout remainder of text 
(xg, Myg, vz n) scaled coordinates 
dimensionless coordinates in units of s 
(x,y>/K, Z% /K) 
sweep angle of leading edge 

Vl -M 2 

D,/D 2 

ratio of specific heats for air 
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At p 
8 


5 1 

5 2 


e 


£,r?,r 


\ 

x i 

P,v 

P 

T 

4> 

v 

x 

* 0 

* 

co 


Jump in \p j at plane of wing or vortex wake 
thickness ratio or measure of camber and angle of attack 

x 2 ' X 1 

*max ^max"^ 

(5M) 2 / 3 

swept wing scaled coordinates (£ = x - y sin ot) 
sin a 

wM/(l -M 2 ) 

scale factors on yQ and zq, respectively (p = v = 5 1 / 3 M 2 / 3 ) 

density made dimensionless to free-stream value 

parameter in analysis of differential operator 

dimensionless perturbation velocity potential 

scaled perturbation velocity potential, 0 = e < p 

0 X + icoi// acceleration or pressure potential 

gicur/r for planar wing , Hq^ 2 \Xt) for airfoil 

^iMXjx ^ fundamental source solution for integral equation 

angular reduced frequency (s x frequency/uQ) 
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DESCRIPTION OF THE METHOD FOR COMPUTING 
UNSTEADY FLOW FOR THE HARMONIC MOTION OF A WING 


A detailed mathematical derivation of the method of finite differences applied to the solution 
of the unsteady velocity potential for the flow about a harmonically oscillating wing is presented in 
the appendix. Formulas are derived for the rectangular wing, for the swept wing, and for 
two-dimensional airfoils. Boundary conditions on the wing, on the trailing vortex sheet, and on the 
outer boundaries of the mesh are treated in sufficient detail to be coded for implementation on an 
electronic digital computer. A brief discussion will be presented here and will be limited to the 
oscillating airfoil, since only the two-dimensional problem was programmed to test the method. 

The complete nonlinear differential equation was simplified by assuming the flow to be a small 
perturbation from a uniform stream near the speed of sound. The resulting equation for steady flow 
is (see appendix, equation 12) 

[K-(t- lV t -(7+ l)^ x ]v? xx + ^ y y-(2^ xt + ^ tt )/e = 0 (1) 

9 9 

where K = (1 - M z )/M z e, M is the free-stream Mach number of velocity Uq in the x-direction, x 
and y are made dimensionless to the semichord s of the airfoil, and the time t to the ratio s/Uq. 
With the airfoil shape as a function of time defined by the relation 

y 0 = 5f(x, t) 


the linearized boundary condition becomes 

Vy — f x ( x ? t) + fj(x, t) (2) 

The quantity 5 is associated with properties of the airfoil such as maximum thickness ratio, 
camber, or maximum angle of attack and is assumed small. The coordinate y is scaled to the 
physical coordinate yQ by the factor 


y = S^M^yQ 


and e is given in terms of 5 by 

e = ~(<S/M) 2 / 3 


8 



The pressure coefficient is found from the relation 

Cp = - 2e (<p x + 1 p t ) 

JThe preceding differential equation is simplified by assuming harmonic motion and by 
assuming the velocity potential to be separable into a steady-state potential and a potential 
representing the unsteady effects. We write for a perturbation velocity potential 

<P = <Pq tx, y) + >p l (x, y)e lcjt 

and for the body shape 

f(x, t) = fg(x) + fj(x)e la,t 

Since the steady-state terms must satisfy the boundary conditions and the differential equation in 
the absence of oscillations, we obtain , 

j 

[K - ( T + 1) v? 0x ] <^> 0xx + = 0 (4) 

with 

^0y = f 0( x )> y = 0,-l <x< 1 (5) 

j 

On the assumption that the oscillations are small and products of may be neglected, equations 
(1) and (2) with the aid of equations (3) and (4) yield 

{[K-(7+ lVo x ]^lx} x + ^lyy-( 2ioj / e )^lx + c l^i =0 (6) 

where 

q= co 2 /e-ico(7-l)<p 0xx 

subject to the wing boundary conditions 

<Ply = f'l( x ) + iwfi(x) = F(x), y = 0, -1 < x < 1 

A computer program for solving the steady-state transonic flow about lifting airfoils based on 
equations (4) and (5) was developed by J. A. Krupp (ref. 4) and is designated as program TEA-330 
by The Boeing Company. The output of this program was used in computing the coefficients for 
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the differential equation of the unsteady potential. The similarity of the unsteady differential 
equation with the steady-state equation suggests that the method of column relaxation used by 
Krupp (ref. 4) for the nonlinear steady-state problem should be 'an effective way to solve equation 
(6) for the unsteady potential <P j. Note that equation (6) is of mixed type, being elliptic or 
hyperbolic whenever equation (4) is elliptic or hyperbolic. Central differencing was used at all 
points for the y derivative and at all subsonic or elliptic points for the x derivatives. Backward (or 
upstream) differences were used for the x derivatives at all hyperbolic points. 

The boundary condition that the pressure be continuous across the wake from the trailing edge 
was found in terms of the jump in potential A<P] to be - 

= A<pj.e - ^ w ( x-x t) 


where A <P^ is the jump in the potential at x = x^. just downstream of the trailing edge and is 
determined to satisfy the Kutta condition that the jump in pressure vanish at the trailing edge. The 
quantity A <Pj is also used in the difference formulation for the derivative ^jyy to satisfy continuity 
of normal flow across the trailing edge wake. 

For the set of difference equations to be determinate, the value of ‘Pj or its derivative must be 
prescribed on the mesh boundary. Following Klunker (ref. 8), we found an asymptotic integral 
representation for the far-field <P j potential and for the related pressure potential <Pj x + iw^j. 
Because of the difficulty with convergence of the integral over the wake for the integral equation of 
the velocity potential, upstream and downstream boundary conditions for the mesh were given in 
terms of the pressure potential ^j x + icu<Pj for which the wake integral can be integrated in closed 
form. The value of was computed at one point on the upper boundary and one point on the 
lower boundary which were conveniently chosen to facilitate rapid convergence of the wake 
integral. The values of </>j at other points on the upper and lower boundaries were found by 
integrating numerically with respect to x the quantity <Pj x + icu<Pj. 

The complete boundary value problem for the oscillating airfoil is illustrated in figure 1 . The 
notation used in the computer program is also shown. Finer mesh in the y direction was used in 
the vicinity of the wing plane. Finer grid in the x direction was used to obtain greater accuracy 
near leading and trailing edges where velocity gradients are higher. 
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i p prescribed from far-field solution 



Figure 1.— Boundary value problem for the unsteady perturbation velocity potential. 



DISCUSSION OF COMPUTED RESULTS FOR THE 
HARMONICALLY OSCILLATING FLAP ON AN AIRFOIL 


The numerical procedure for computing the transonic unsteady flow described in the foregoing 
section was applied to the airfoil with the harmonically oscillating flap. For the basic steady 
transonic flow, we chose the flow about a NACA 64A006 airfoil at zero angle of attack for 
free-stream Mach numbers of M = 0.8 and 0.85. The steady flow potential <Pq used in computing the 
coefficients for the unsteady method was computed by the program TEA-330 developed by 
J. A. Krupp (ref. 4). The distribution of pressure coefficients on the airfoil for the two Mach 
numbers is shown in figures 2 and 3. Mach numbers and the airfoil were chosen for the purpose of 
comparing computed results with the experimental measurements of Tijdeman and Schippers in 
reference 9 for an oscillating flap of quarter-chord length. We note that for M = 0.8 the airfoil is 
subcritical, while for M = 0.85 it is supercritical over a small region. The reduced frequencies for the 
flap oscillating at 30 cycles per second in free-stream Mach numbers of M = 0.8 and 0.85 for the 
wind tunnel speeds used by Tijdeman and Bergh in reference 9 are w = 0.060 and 0.056, 
respectively. The grid consisted of 75 points along y = constant lines and 58 points along 
x = constant lines, yielding a total of 4350 points in the mesh. The far-field boundaries are at 
x = -2.62 and 2.75 and y=±6. This same grid was used for both the steady and unsteady 
calculations. 

The unsteady contribution to the location of the airfoil surface for the harmonically oscillating 
flap as a function of chord position is given by 

\ 

f j (x) = 0, -1 < x < 1 - L f 

f 1 (x) = a(x- 1 + L f ), 1 - L f < x < 1 

where the airfoil is defined by -1 < x < 1 and Lf is the length of the flap in units of semichord. For 
the experiments of Tijdeman and Schippers (ref. 9), Lf = 1/2, and the amplitude was 1.5° The 
boundary conditions in equation (82) of the appendix become 


•Ply - 0> -1 < X < 1 - Lf, y - 0 ( 7 ) 

^ly = (a/5) [1 +icu (x - 1 + Lf)] , 1 - Lf <x < 1 
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Figure 2.— Distribution of coefficient of pressure, C^, at Mach 0. 794 



Figure 3.— Distribution of coefficient of pressure, C^, at Mach 0.853. 
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Thus, the quantities and F-^\ defined at the mesh points Xj and used in the difference 

equations near the boundary, become 


F XU) = F .(L) = (a/g)[ ! + i0J (x . _ i + Lf )] , i > i a , 

F .(U) = F .(L) = 0 , i<ia , 


where i a is the x-index in the mesh satisfying the relation 


< 1 -L f <x 


*a+l 


' To check out the program, the flow at Mach number M = 0.8 over a flat plate was computed 
for the oscillating boundary conditions of equation (7) and for u = 0.06. For this case 
<Pq x = ^Oxx ~ 0- The resulting differential equation is equivalent to the classical linearized unsteady 
subsonic differential flow equation for harmonic motion; namely, 

(1 - M^)<^ xx + <£>yy - 2iM^c o<^ x + = 0 

This equation has been studied extensively and several kernel function integral equation methods 
have been developed to solve for the related acceleration or pressure potential, <Pj x + icu^j. 

The jump in pressure coefficient computed by the present method for the flat plate at M = 0.8 
and u) = 0.06 is given in figures 4 and 5 for the real and imaginary parts, respectively. The 
convergence of the solution was very slow. Overrelaxation with a parameter of 1.3 in the manner of 
Murman and Cole (ref. 2) and Krupp (refs. 3 and 4) was used in the iteration process, and the 
relaxation calculations were made by sweeping alternately downstream and upstream for many of 
the iterations. No definite advantage in alternate direction sweeping of the field as compared with 
the customary mode of sweeping downstream followed in references 2 through 4 was proven, 
although Newman and Klunker (ref. 6) found that symmetry of nonlifting solutions for steady 
subsonic flow was preserved. For the far-field boundary conditions, the value of the unsteady 
potential <P j was computed for upper and lower mesh boundaries, and < Pj x + iaj'^j boundary 
conditions were prescribed following the discussion in the section on far-field boundary conditions. 
The far-field boundary data were updated every 100 iterations near the beginning of the 
calculations and every 20 iterations for the later calculations. The more frequent updating seems to 
be more efficient in obtaining the final result. With all values of <Pj in the mesh starting with zero, 
the complete calculations shown in figures 4 and 5 were accomplished in about 3000 iterations to a 
convergence of 2 x 1 0'~> for the magnitude of the maximum difference between iterations. In spite 
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Figure 4.— Real component of the jump in unsteady pressure coefficient across a flat plate 
with a harmonically oscillating quarter-chord flap. Free-s tream Mach number 
M = 0.8 and reduced frequency co = 0.06. 




- Present method 

- Collocation method of Rowe (ref. 10) 



x 


Figure 5.— Imaginary component of the jump in unsteady pressure coefficient across a 
flat plate with a harmonically oscillating quarter-chord flap. 

Free-s tream Mach number M = 0.8 and reduced frequency c o = 0. 06. 

of the complex arithmetic and the more complicated equations, the computing time for each 
iteration is about 1 second on the CDC 6600 compared with the steady flow program which 
requires about 2/3 second per iteration for the 58 by 75 grid. 

The results are compared in figures 4 and 5 with calculations by the kernel function of 
reference 10 for three-dimensional flow at the center of a rectangular wing of aspect ratio 20. The 
real parts of the jump in pressure coefficient by the two methods agree quite well, but the 
imaginary part from the present method gives amplitudes with less negative value near the leading 
edge and larger positive value on the flap. The x,y mesh points were those used to compute the 
steady-state potential using TEA-330 and contain a finer spacing near the leading and trailing edges. 
The grid points used in the calculations are indicateu by circles on graphs of computed results. It 
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has been shown in reference 1 1 that there is a logarithmic singularity at the hinge point and an 
infinite square root singularity at the leading edge. Because of these singularities, no mesh points for 
the unsteady flow should lie on the lines x = -1 and 0.5. Increasing the fineness of the x-grid in the 
vicinity of x = 0.5 had no noticeable effect on the solution for the^flat plate, however. The results 
also appear to be insensitive to the far-field boundary conditions. Setting <Pj x + iu></>j =0 on all 
four far mesh boundaries located at y =±6.5 and x =±3.5 yielded differences scarcely observable in 
figures 4 and 5. 

The solution for the flat plate with a free-stream Mach number M = 0.8 and reduced frequency 
U) = 0.06 was used as a starting solution for u) = 0.1794, corresponding to a frequency of 90 cycles 
per second. The singularity at the leading edges for u) = 0.1794 has the opposite sign from that at 
co = 0.06. After 2800 iterations the solution converged to 10"^, and the real part along with the 
equivalent solution from the collocation method of Rowe is shown in figure 6. The agreement is 
quite good. Similarly, the comparison of the finite difference solution for the imaginary part of the 
pressure jump with the corresponding solution by the collocation 'method is shown in figure 7. The 
agreement is better than that for u> = 0.06 (fig. 5). The rate of convergence is illustrated in figures 8 
and 9. The cause of the large number of iterations required for convergence is apparently the change 
in sign of the leading edge singularity for co changing from 0.06 to 0.1794. 

With the calculations of the flat plate solution for M = 0.8 and cv = 0.06 as starting values for 
<p j, the solution of the complete linear differential equations for Mach 0.794 and tu = 0.064 with 

u = K - (7 + l)^0x • 

provided by the solution from TEA-330 was obtained in about 1200 iterations. Agreement with the 
experimental results of Tijdeman and Schippers (ref. 9) is not as good as the agreement between the 
classical linear solution and the present method for flat plate flow. The pattern of the present 
method agrees with that obtained by Tijdeman and Schippers (ref. 9), as shown in figures 10 and 
11. The computed values of real ACp are considerably larger than the experimental measurements. 
The discrepancy between experiment and theory is not understood. There are a 
number of possibilities: boundary layer or separation effects, or both, near 

the aft portion of the airfoil, or unknown problems associated with the theory 
or with the pressure measurements. However, the validity of the calibration of 
the pressure tube system used by NLR is confirmed by Destuynder and Tijdeman in 
reference 12, and should provide reasonable experimental results. 

Starting with the M = 0.8 and U) = 0.06 mesh values of <Pj and with^Q from TEA-330, we 
computed the solution for M = 0.853 and u; = 0.060 with another 1200 iterations to a maximum 
error between iterations of about 3.5 x 10‘^. The results for the jump in pressure coefficient on the 
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■g oo Present method 

Collocation method of Rowe (ref. 10) 



Figure 6.— Real component of the jump in unsteady pressure coefficient across a flat 
plate with a harmonically oscillating quarter-chord flap. Free-stream ! Mach 
number M = 0.8 and reduced frequency c o = 0 1 794. 




X 


Figure 7.— Imaginary component of the jump in unsteady pressure coefficient across a 
flat plate with a harmonically oscillating quarter-chord flap. 

Free-stream Mach number M = 0.8 and reduced frequency co = 0.1794. 


airfoil are shown in figures 12 and 13 for the real and imaginary parts, respectively. The pressure 
peak for both the real and imaginary parts close to the midchord as observed by Tijdeman and 
Schippers (ref. 9) is also predicted by the present method. Although magnitudes are somewhat 
higher, the pattern of the theoretical imaginary part from the present method agrees with the 
experiments. 


By comparing the two solutions in figures 10 through 13, one can see the influence of 
increasing Mach number on the propagation of disturbance from the oscillating flap. Note that the 
real component of the pressure amplitude decreases with Mach number near the leading edge. A 
peak builds up at the semichord point due to the local supersonic region, which becomes subsonic 
just downstream of the midchord. These observed properties are common to both the theoretical 
and experimental results. 

The solution for cj = 0.06 and M = 0.794 was used as the starting solution for w = 0.253 and 
M = 0.804, corresponding to a frequency of 120 cycles per second. After 4000 iterations the 
solution was converged to 4 x 1 0"^ for the maximum difference between iterations. The column 
relaxation was swept alternately upstream and downstream and an overrelaxation parameter of 1 .4 
was used. This flow is completely subsonic. The pattern of the real part of the jump in pressure 
coefficient from the finite difference method in figure 14 matches the experimental data as in the 
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Figure 8. —Rate of convergence of the real part of the jump in pressure coefficient on a 

flat plate with oscillating quarter-chord flap for a change in reduced frequency 
from co = 0.06 to lo = 0.1 794 and free-stream Mach number M = 0.08. 


other cases, but here the magnitude is in somewhat better agreement. The imaginary part shown in 
figure 15, however, has much higher negative values near the leading edge than the experimental 
data, but the match is better aft of the hinge line. 

The solution for u> = 0.056 and M = 0.85 was also extended to oj = 0.179. This solution was 
run 3600 iterations and the convergence is 1.93 x 10"^ maximum difference between iterations. 
From figure 16 just downstream of midchord, a peak in the real part of the jump in pressure 
coefficient has developed and increased to 12, while the imaginary part in figure 17 shows a peak of 
-20 just upstream of midchord. The experimental measurements also show these peaks, but on a 
much smaller scale. The amplitude buildup occurs near the sonic line in the decelerating region of 
the flow for both values of reduced frequency (see figs. 12 and 13). From the experience in 


20 










Present method 

Experimental results of Tijdeman 
and Schippers (ref. 9) 



Figure 10. —Real component of the jump in unsteady pressure coefficient across 

NACA 64A006 airfoil with a harmonically oscillating quarter-chord flap. 
Free-stream Mach number M = 0.794 and reduced frequency go = 0. 064. 




Figure 1 1.— Imaginary component of the jump in unsteady pressure coefficient across 
NACA 64A006 airfoil with a harmonically oscillating quarter-chord flap. 
Frees tream Mach number M = 0. 794 and reduced frequency co - 0. 064. 
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Figure 1 3.— Imaginary component of the jump in unsteady pressure coefficient across NACA 64A006 airfoil with a harmonically 
oscillating quarter-chord flap. Free-stream Mach number M = 0.853 and reduced frequency co = 0.060. 




Figure 14.— Real component of the jump in unsteady pressure coefficient across 

NACA 64A006 airfoil with a harmonically oscillating quarter-chord flap. 
Free-stream Mach number M = 0.804 and reduced frequency co = 0.253. 





Figure 15.— Imaginary component of the jump in unsteady pressure coefficient across 
NA CA 64 A 006 airfoil with a harmonically oscillating quarter-chord flap. 
Free-s tream Mach number M = 0. 804 and reduced frequency oj = 0. 253. 





Figure 1 6. —Real component of the jump in unsteady pressure coefficient across NACA 

64 A 006 with a harmonically oscillating quarter-chord flap. Free-s tream Mach 
number M = 0. 858 and reduced frequency oj = 0.1 79. 
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Present method 

Experimental results of Tijdeman 
and Schippers (ref. 9) 


Figure 1 7.— Imaginary component of the jump in unsteady pressure coefficient 

across NACA 64A006 airfoil with a harmonically oscillating quarter-chord 
flap. Free-s tream Mach number M= 0. 858 and reduced frequency co = 0.1 79. 




calculation with the amplitude of oscillation of 1-1/2°, it appears that calculations must be 
continued until the maximum magnitude of the difference between iterations is about 1CT'\ 

From a study of figures 1 0 through 1 7, we see that solutions from the finite difference method 
predict the essential features of the pressure distributions found by experiment, whereas the 
linearized subsonic kernel function method does not. The larger amplitudes may be due to 
inaccuracies in the method of measurement, to boundary layer and separation effects on the aft 
portion of the airfoil, or to wind tunnel wall influence. 


RECOMMENDATIONS FOR INCREASED ACCURACY AND REDUCED COMPUTING TIME 


It is pointed out in the appendix that the difference equations for the flow field result in a 
large set of linear equations for the values of the unsteady velocity potential at the mesh points. 
Column relaxation methods were found to require a large number of iterations to obtain the 
solution. Since the matrix of the equations contains only five or six nonzero diagonals for 
two-dimensional flow (and only seven or eight for planar wings), the matrix is very sparse, 
containing mostly zeros. 

Special methods which take advantage of this property should be applied to solving this set of 
equations. One possible technique is the Gram-Schmidt orthogonalization method of references 1 3 
and 14, which construct a solution from a set of vectors orthogonal to the set of vectors spanned by 
the coefficients and right-hand side terms of equations (127) and (128). Scalar products used in 
constructing these orthogonal vectors wil) contain only six or seven products at most. The actual 
number of operations should be considerably less than the 1 200 to 3000 iterations of the column 
relaxation required for the convergence by the column relaxation method. As formulated in this 
manner, the complete calculations will have to be carried out several times to match far-field 
boundary conditions and to satisfy the Kutta condition. Each updating of the far-field boundary 
conditions and the Kutta conditions from the last solution requires only a modification of the 
right-hand side terms in preparation for the next solution. Satisfying the Kutta condition could be 
incorporated into the solution of the linear equations, but this would involve additional nonzero 
elements in the matrix for those points adjacent to the plane of the wake. Many terms in the matrix 
are duplicated, and taking advantage of this feature can reduce some of the computation. A study 
of the various special methods available for solving large orders of linear equations should be made 
to determine which method is best suited to this particular problem. 
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It is also possible to improve the efficiency and accuracy of the column relaxation procedure. 
This may take the form of finding optimum values of the overrelaxation parameter at subsonic 
points and of underrelaxation at supersonic points, or possibly changing the mesh size near the 
leading edge and near the hinge line. 

In calculating the exact potential steady transonic flow over the blunt airfoil, Garabedian and 
Korn in reference 1 found that obtaining first a solution with a coarse grid and then refining the 
grid to obtain greater accuracy reduced significantly the number of iterations required for 
satisfactory accuracy in the final solution. This procedure should improve both efficiency and 
accuracy of the unsteady finite difference solution as well. 


Experience in computing additional unsteady flows, such as harmonic pitching of the airfoil 
and comparison with careful measurements on an oscillating airfoil, is needed to clarify the validity 
of the method. 


EXTENSION OF THE WORK TO THREE-DIMENSIONAL WINGS 


In the appendix, a complete analysis based on transonic small perturbation theory is given for 
applying the method of finite differences to the solution of unsteady flow about finite wings of 
both rectangular and swept planforms. All formulas required for programming the method are 
presented. Much of the computer logic in the airfoil program could be directly employed in the 
program for finite wings. The basic approach would be essentially unchanged, but with the 
additional dimension adding more complication and a greater number of points. The program either 
of Newman and Klunker (ref. 6) or Ballhaus and Bailey (ref. 7) would be required to provide the 
values of the steady perturbation potential at the grid points for the coefficients of the unsteady 
flow difference equations. 


EXTENSION OF THE WORK TO ELIMINATE THE SMALL DISTURBANCE ASSUMPTION 


' ' '• The general procedure followed here for setting- up the linear equation for the unsteady 
velocity potential could be applied to the more complete potential equation. Garabedian and Korn 
(ref. 1 ) transformed the independent variables of the differential equation by a conformal mapping 
of the physical plane to map the flow field about a blunt airfoil into the interior of a unit circle. 
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The transformed differential equation is expressed in difference form in the new variables and 
solved numerically by relaxation in the manner of Murman and Cole (ref. 2) and Krupp (refs. 3 and 
4). In this way the airfoil boundary conditions are satisfied exactly and far-field boundary 
conditions are easily prescribed. 

This same transformation of the independent variable could be similarly applied to the 
differential equation for the velocity potential of unsteady two-dimensional flow. By representing 
the complete flow as the sum of a steady-state potential and an unsteady potential, the transformed 
differential equation for the unsteady flow could be linearized in the same fashion by using the 
condition that the steady-state potential satisfies the nonlinear equation and by neglecting all terms 
containing products of the unsteady potential and its derivatives. The unsteady boundary 
conditions could be linearized by satisfying them on the boundaries of the stationary airfoil and 
steady wake. The techniques for solution would differ more in the increase of actual computation 
required for the coefficients than in the actual number of difference equations to be solved. 

The Boeing Company has developed a similar method for solving the complete nonlinear 
differential equation for the steady transonic flow potential about a blunt airfoil by using as 
independent variables the velocity potential and stream function from the incompressible solution. 
With the differential equation for steady flow expressed in these variables, the boundary conditions 
on the airfoil are satisfied exactly. The same procedure could also be applied to this formulation. 
There appears to be no basic reason why the numerical procedure used for the differential equation 
based on the assumption of small perturbations near sonic velocity could not be carried over to the 
blunt airfoil solutions. 

Since no mapping technique exists for simplifying the application of the boundary conditions 
on three-dimensional wings, the problem of eliminating the restriction to small perturbations for 
transonic flow is considerably more difficult. A fairly elaborate mesh system must be constructed 
near the wing surface, and a method of satisfying the boundary conditions must be devised. 


CONCLUSION 


A method for applying the finite difference technique to the solution of the unsteady 
transonic small perturbation differential equation has been developed analytically with sufficient 
detail in this document that it can be programmed for computation on a high-speed digital 
computer. Coefficients of the differential equation contain the nonlinear steady flow solution about 
the wing and must be computed and supplied to the unsteady flow program at the outset of 
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computations of the unsteady flow. To test the method, the technique for the oscillating airfoil was 
programmed and computations were made for the oscillating flap of quarter-chord length on a 
NACA 64A006 airfoil. 

Analysis of the truncation terms of the difference equation indicates that the stability of the 
difference formulation for subsonic points is not dependent upon the relative grid sizes in the three 
directions. For supersonic points, a viscous-like term was found in the truncation terms which is 
similar to that found by Murman and Cole (ref. 2) to yield stable and smooth calculations near 
shock waves. Calculations using the method indicate that the procedure is indeed stable and 
convergent. 

Five cases were computed involving two free-stream Mach numbers and two frequencies. The 
numerical method was found to be stable for flows with local supersonic regions as well as strictly 
subcritical flows. The pattern of the distribution of pressure coefficient on the airfoil in all cases 
follows closely the distributions observed by Tijdeman and Schippers (ref. 9); however, the pressure 
levels in the calculated cases are generally, higher than those found by Tijdeman and Schippers (ref. 9). 

The discrepancy between theory and experiment may be due to boundary layer and 
separation effects near the aft portion of the airfoil, possible wind tunnel 
wall interference, or to unknown problems associated either with the theory or 
with the pressure measurements. The validity of the calibrations used by Tijdeman 
and Schippers (ref. ,9) is confirmed by Destuynder and Tijdeman in reference 12. 

For 75 grid points in the x or stream direction and 56 in the y direction, the number of 
iterations for convergence to a maximum difference between iterations of 2 x 10"^ varied from 
1200 to 3000 depending upon the starting solution. Each iteration takes about 1 second on the 
CDC 6600 as compared to 2/3 second for the steady-state program TEA-330 developed by Krupp 
(ref. 4). More computations are needed to improve the rate of convergence by finding optimum 
values of the overrelaxation parameter for subsonic points and of the underrelaxation parameter for 
supersonic points, and to establish suitable criteria for acceptable convergence of the solution. Since 
Garabedian and Korn (ref. 1) found that first obtaining a solution with a coarse grid and refining 
the grid to carry on the calculations reduced significantly the number of iterations required for 
solution of the steady-state equation, a similar procedure should also reduce computing time for the 
unsteady problem. See note added in proof on page 34. 

Efficiency also may be improved by solving the complete set of equations instead of the 
column relaxation used in the study. Because the set of difference equations is linear in the 
unsteady potential and because the matrix of coefficients contains only a few nonzero diagonals, a 
special procedure for solving the complete set of equations at once seems feasible and should be a 
subject of future study. 
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The present method was formulated on the assumption of small perturbations near the speed 
of sound, but the method could be extended without too much difficulty to tfce flow around a 
blunt airfoil for which the steady flow boundary conditions are not linearized but are satisfied 
exactly. The two-dimensional problem is made easy by the possibility of using body coordinates to 
simplify the application of boundary conditions in the finite difference method. Since no simple 
body coordinate system can be formulated for the three-dimensional wing, a more elaborate mesh 
system would be required in the neighborhood of the wing surface, making it more difficult to 
eliminate the restriction of small perturbations. The unsteady boundary conditions presumably 
would be linearized and applied at the steady-state wing position. 

The results of the present study indicate the feasibility of the proposed method of finite 
differences for finding the unsteady transonic flow about harmonically oscillating wings. 


Boeing Commercial Airplane Company 
P.O. Box 3707 

Seattle, Washington 98104, February 25, 1974 


Note added in proof: 

Recent studies have shown that, by the use of larger values of the over- 
relaxation parameters such as 1.8 to 1.9, the number of iterations to obtain a 
converged solution was reduced by an order of magnitude when applied to a re- 
computation of one of the flat plate examples. 


34 



APPENDIX 


Derivation of the Differential Equation and Boundary 
Conditions for Harmonic Motion of the Thin Wing in Transonic Flow 


Let 


u = u 0 ( 1 + 0 x o >» v = u O 0 y o > w = u 0 % 


( 1 ) 


be the velocity components of the flow, where 0 is the perturbation velocity potential, XQ,yQ,zg 
Cartesian coordinates, tg the time, and ug the undisturbed flow velocity. As given in Landahl 
(ref. 15), the equations of motion are 


% + [p(l +0 XQ 


)]xn + ^^yg ) y O + (p 0 Zn ) z r = 0 ’ 


z o' z o 


( 2 ) 



1/M 2 - ( T - l)[0 t + 0 X() + 



(3) 


where p is the density, c the velocity of sound, 7 the ratio of specific heats, and M the free- 
stream Mach number. The coordinates (xQ,yQ,zg) are made dimensionless to the wing semichord 
and the time tg to the ratio of the semichord to the velocity ug. Since 

dc 2 /(7- 1 ) = c 2 dp/p , ( 4 ) 


we have from equation (3), 

c2p t 0 /p = ‘ [ 0 tgtg + 0 Xgtg + ^Xg^xgtg + fy/ygtg + ^Zg^Zgtg ] 


(5) 


Similar relations for c 2 P X g/p, c ^Pyg/P> and c 2 P Z g/p may be written down by replacing the last 
subscripted variable tg with xg, yg, and zg, respectively. By expanding the derivatives in equation 
(2), eliminating c 2 Ptg/P by equation (5), and eliminating c 2 p X q/p, c ^Pyjp, and c 2 P Z g/p by 
equations similar to equation (5), the following differential equation for 0 results 
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( 6 ) 


Vo 4 " **o*o + V*o*o + 0 yoVo + 

+ [ 0 t o x o + 0 x o x o + 0 x o 0 x o x o + 0 yo 0 yo x o + 0 z o 0 z o x o] C + 0 x o) 

+ [ 0 t oyo + 0 x oyo + 0 x o 0 x c>yo + 0 yo 0 yoyo + 0 z o 0 z oyo] 0 yo 
+ [ 0 t o z o + 0 x o z o + 0 x o 0 x o z o + 0 yo 0 yo z o + ' 0 z o 0 z o z o] 0z o 
' [i/M 2 • (T ■ 1) ( 0 t o + 0 X O + 0 x o 2/2 + 0 y o 2/2 + 0 z o 2/2 )] ( 0 X O X O + 0 yoyo + 0 z o z o) = 0 

Since we shall consider small perturbations from a parallel flow, it is convenient at this stage to 
discard third-order terms in equation ( 6 ). Combining terms and simplifying yields 

[l - M 2 - (7 ■ 1 )M 2 0 tQ - (7 + 1 )M 2 0 Xq ] 0 x o x o + 0y o y o + 0 Z O Z O 
- 2M ^x 0 t 0 - M 2 0 t() t o - m2 [ 2 ^x/x 0 t 0 + 2 Vyoto 

(7) 

+ 20 z o 0 z o t o + 2 , % 0 yo x o + 20 z o 0 z o x o 
-(T- 1 > (% + % ) ( 0 y o y o + 0 z o z o) ] = 0 


We now introduce a small parameter e and scale factors ju and v on yg,zg and, by an order 
of magnitude comparison of terms of equation (7), we determine which second-order terms in 
equation (7) are to be retained for free-stream Mach numbers near 1 . Let <p = e< p, x = Xq, t = tg, 
y = /uyg, and z = vzq. The quantity e is a small quantity whose magnitude in terms of the thickness 
ratio is to be determined later in the development. The quantities ju and v are scale factors for 
setting the orders of magnitude of the principal linear terms. With these substitutions the differ- 
ential equation becomes 

e[(l - M 2 ) - (7 - 1 )M 2 e^ t - (7 + 1)M 2 e^ x ] ip xx + p 2 e 0 yy 
+ t’ 2 e<^ zz - 2M 2 e<^ xt - M 2 e<^ tt - M 2 e 2 [2v? x v? xt 
+ 2M 2 vy> yt + 2v \z^t + + 2v2 *z*zx 

+ (7 - l)(¥>t + ^x^^yy^ 2 + ^zz ^ = °- 
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Since our primary interest is in transonic flow, we assume that the small quantity e is of the same 
order as 1 - M 2 . Dividing the equation by M 2 e 2 and setting (1 - M 2 )/M 2 e = K leads to 

[K - (t - 1 V t - (7 + 1 V x ] ^xx + (M 2 /M 2 eV yy + (v 2 /M 2 e>M zz 
- (2* xt + * tt )/e - [ 2^ xt + 2/x 2 v? y ^ yt + 2v 2 * z ^ t (8) 

+ 2n 2 ip y <p yx + 2v 2 V z <P zx + (7 - 1)0P X + ^t^ 2| fyy +i;2 ^zz )] = °- 

It is convenient to set v 2 /M 2 e = 1 and p 2 /M 2 e = 1 , thereby making y? yy and <p zz terms of the same 
order as the <m xx term. Equation (8) then becomes 

[K - (7- 1)^ - (7 + 1 V x ^xx + {fi yy +l ^zz 

-(2<p xt +« tt )/e-^2<p x v xt + 2v 2 [v y v yt + v z y> zt ( 9) 

+ ^y^yt + ^z^zx + ( 7- 1 XsPt + ^x^^yy + V zz )/ 2 ] | = 0. 


To determine m in terms of the thickness ratio 5 we consider the boundary conditions on the 
wing; we let 

F(x 0 , y Q , t Q ) = z 0 - 5f(x 0 , y Q , t Q ) = 0 


describe the wing surface at each time tQ. Then the boundary conditions that the flow be tangential 
to the wing surface become 

DF/Dt=0 


or 




Y ) * 5f. 0 + 0 

x o yo yo 2 


5f t _ = 0. 


0 


( 10 ) 


In scaled coordinates with 0 = e<p, we have 

* z - (5/e/i)(f x + f t ) - (S/ M )f x * x - (5ju )f yV3y = 0. 
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Letting 5/eju = 1 makes all the linear terms have the same order of magnitude. Then the boundary 
conditions become 


' f X * f t - 6 ^X f X ' ^ W = 0. 


(ID 


9 9 

Since p = 5/e and )U z /M z e = 1, we can solve for /a and e in terms of 5 . Thus 


and 


e = (5/M) 2 / 3 


= v = 5 !/ 3 M 2 / 3 = Me ! / 2 


Upon neglecting terms of the order of e in equation (9), we finally obtain the differential equation 


[K-(t- lV t *(7+ l)^ x ^xx + ^yy + ^zz 

- 2y> x y> xt - (2<p xt + v? tt )/e - 0. 


( 12 ) 


The boundary condition from equation (11) 

^z = f x + f t 

is to be satisfied on the wing, represented by its projection on the zq = 0 plane. 

For convenience, we separate the perturbation potential <p into a steady flow ipg and an 
unsteady flow ipj. Thus, we let 


V = <Po( x ’ y ’ z ) + ¥>i(x,y,z,t). (13) 

The steady flow term y?Q must satisfy equation (12) independently. Thus, we obtain 

[K-(t+ lV 0x ]^0xx + ^0yy + ^0zz = 0 - (14) 

This is the well-known transonic small perturbation equation for steady flow studied recently by 
Murman and Cole, Krupp, Newman and Klunker, and Ballhaus and Bailey (see refs. 2 through 7). 

It is nonlinear and of mixed type. When <pg x < K /(7 + 1), the equation is elliptic; and when 
y>Q x > K /(7 + 1), it is hyperbolic. The quantity K /(7 + 1) is seen to be the sonic velocity for the 
flow represented by the solution of equation (14). 
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Substituting equation (13) into (12), making use of equation (14), and neglecting product 
terms in yields the following linear differential equation for the unsteady perturbation potential: 


[K-(7+ D^ 0x ]^l xx - [(7- 1 )^it + (T+ D^lx^Oxx 

+ ^lyy + ^lzz - ^lxt + ^ltt)/ e = °- 


(15) 


Here we have also neglected the coefficient i£>q x of the i^j xt term as negligible compared to 1 /e. 
We shall assume harmonic motion for the wing; then with 


*\ ~ ¥>i(x, y, z)e 


icot 


(16) 


equation (15) becomes 

{ [K-( 7 + l)^Ox^lx}x + ^lyy + ^lzz'( 2 ic °/ e ^lx + '(?* Di^Ozz^l = 0 
This may be written as 


[u^ lx -(2iw/eV 1 ] x + ^. lyy + ^ lzz +qv» 1 = 0 (17) 

where u = K - (7 - 1)<Pq x and q = co 2 /e - icc (7 - l)<Poxx - Equation (17) is linear and like equation 
(14) is of mixed type, being elliptic or hyperbolic whenever the solution of equation (14) is elliptic 
or hyperbolic. 


Formulation of the Difference Equations for Numerical 
Solution of Unsteady Flow From Harmonic Motion of a Thin Wing 

Since tpQ is a numerical solution of equation (14) whose values are known only at discrete points 
of a mesh, equation (17) is also to be solved numerically. For this purpose equation (17) is expressed 
in difference form. For subsonic flow points, we use central differences for derivatives in the x,y,z 
directions. For hyperbolic or supersonic points we use backward (or upstream differences) in x but 
retain central differences in the y and z directions. This is the practice followed by Murman and 
Cole (ref. 2), Krupp (refs. 3 and 4), Steger and Lomax (ref. 5), and Newman and Klunker (ref. 6 ), 
since only upstream data can influence a hyperbolic point in regions of local supersonic flow. Thus, 
for subsonic flow at the point Xj, yj, z^, we have 

(u^l x )x - { u i+ V 2 jk^i+ljk "^ijk^ x i+l * X P ' ^-'Ajk^ijk '^i-ljk^ x i ‘ x i-l^} t 2 /( x i+l x i-i ) 1 
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where 


u i+'/2jk - K - (7 + DC^Oi+ljk -^OijkV( x i+l * x i> 

u i-y 2 jk = K - (7 + 1 Xv’oijk - ^Oi-ljk>/( x i - x i-l) 


and i ± l A denotes the value of the quantity at 


x = (x i + x i±1 )/2, 


the x midpoint of the mesh. 

The x derivative may be written more simply as 

( u< ^lx)x = ^ c i u i+ 1 /2jk^i+ljk ' ^ijk^ 

' 2d i u i- 1 / 2 jk (^ijk ' ^i-ljk> 

where 

c i = 1 /( x i+l - x i-l>( x i+l - x i) 
di = 1 / ( x i+ 1 * x i-i X x i • x i-i) 

For the first derivative, we use a second-order difference equation. Writing 

^lx = c lM+ljk*^ijk) + d lMjk’^i-ljk) (20) 

and expanding the terms about the point ij,k, we find the values of c^ and dj j by setting the 
coefficient of <£] x equal to 1 and of y>j xx equal to zero on the right-hand side of equation (20). 
This leads to 


c li = ( x r x i-l) c i 

d li = ( x i+l - x i> d i 
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With these definitions, the complete x derivative term becomes 
I u *l x ’ C2ico/eVi ] x = - 2[ Ci u i+ i /2jk + d iUi .i /2jk 

+ i(cj/e)(d ir c li )] 1 p i j k 

+ 2[cjUj + i^j k - icocjj/e]<£j + jj k 

• + 2[di u i + i /2 jk + icodii/e] v?i_ijk 

Similarly, we write for the central differences of V\yy and <p] zz terms 


or 


where 


2 ^ij+lk ~^ijk ^ijk ~^ij-lk 

^yy'yj+i-yj-iL y j+1 -yj ' yj-yj., 

^lyy = 2a yj^ij-lk ‘ 2 ( a yj + b yP^ijk + 2 ^yj^ij+lk ^ 
*^lzz = 2a zk^ijk-l ' 2 ^ a zk + b zk^ijk + 2E zk^ijk+l 


a yj = 1/(y j+l - y j-l )(y j- y j-D 
b yj = 1/(y j+l ' y j-l )(y j+l - y j>- 


( 21 ) 


( 22 ) 


(23) 


The formulas for a zk and b zk are similar to equation (23) with z and k replacing y and j , 
respectively. Substituting equations (20), (21), and (22) into the differential equation for (pj yields 

-2( E l + E 2>^ijk + 2E l^i+ljk + 2E 2^i-ljk + ^ijk^ijk 

+ 2a yj^ij-lk ' 2 ( a yj + b yjVijk + 2b yj^ij+lk 


+2a zk^ijk-l ' 2 ( a zk + b zkMjk + 2b zky ? ijk+l 


where 


E 1 = c i u i+&jk ' iwc 1 i / e > E 2 = d i u i+ 1 / 2 jk + iwd 1 i/ e - 
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In solving the equations, we relax a column at a time. At each planar point Xj, yj, we solve for 
^ijk’ k = 2 to k max - 1 . Therefore, we place all other terms on the right-hand side and obtain for the 
elliptic point form of the difference equation: 


a zk^ijk-l ‘ ( a yj + b yj + a zk + b zk + E 1 + E 2 - <4jk/ 2 Mjk 
+ b zk^ijk+l = " E Fi+ljk ‘ E 2^i-ljk ' a yj^ij-lk ' b yj^ij+lk 


(24) 


For fixed i j, and for k ranging from 2 to k max - 1 , this is a set of k max - 2 equations to be solved 
for an equal number of <£y k provided that ^ j and ^ijkmax are P rescr ibed by the far-field boundary 
conditions. This yields a system of equations with a tridiagonal matrix for the coefficients which is 
simple to solve. 

When the point is hyperbolic, backwards differencing in x is used. This may be achieved by 
shifting all the i indices one less for the term (u^i x ) x - For yq x we derive a three-point backwards 
difference formula which is second-order accurate in the mesh spacing. For the point i j,k, we write 


*lx “ c 2Mjk ' ^i-ljk) + d 2i^i-ljk * *i-2jk>- 


(25) 


We then determine C2j and d2j so that, when and are expanded in a Taylor’s series in 

x about the coefficient of on the right-hand side equals unity and the coefficient of ^] xx 
vanishes. This leads to 


c 2i = l/(xj - Xj_j) + l/(xj - Xj_ 2 ) 
d 2i = -d 1 i- 1 ■ 


(26) 


The x derivative terms for hyperbolic points then take the form 

[u'Pjx - (2ico/eVj ] x = 2 c m u i .i/ 2 j k (v? i j k - ^-ij k ) 

- 2d i-l u i-3/2jk^i-ljk-^i-2jk) 

- (2ico/e )[c 2i («pyk ' ^i-ljk> ‘ d li-l^i-ljk * ^i-2jk^ • 
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The complete difference equation equivalent to equation (17) in tridiagonal form for hyperbolic 
points becomes 


where 


Wijk-l ' ( a yj + b yj + a zk + b zk " E 3 + c lijk/ 2 Mjk + b zk+l^ijk+l 
= (E 3 + E^j.ijk - E4^i-,2jk|ljk * a yj^ij-lk' b yjfij+lk 


(27) 


E 3 = c i-l u P/2jk- iwc 2i/ e 
E 4 = d i- 1 u i-3 /2jk " iojd 1 i- 1 / e 


Equations (24) and (27) are the difference equations for interior points in the flow field. Since 
boundary conditions on the wing and on the vortex sheet are expressed in terms of the derivative of 
ifil with respect to z , it is convenient to locate the linearized boundary of the wing and wake on 
the plane between k = k m and k = k m + 1 . For these values of k , equations (24) and (27) must be 
modified to incorporate these boundary conditions, and the appropriate relations are developed in 
subsequent sections. 

Boundary Conditions on the Wing 

For harmonic motion, the linearized boundary condition following equation (12) becomes 

^iz =f ix( x >y) + iwf i< x ’y) ( 2g ) 


where we have substituted 


zq/5 = f(x,y,t) = f Q (x,y) + f 1 (x,y)e 10jt 

for the shape of the wing surface with harmonic oscillations. The steady flow perturbation potential 
satisfies the boundary condition i£q z = fg x where zq = 5fg(x,y) defines the undisturbed wing shape. 
The boundary conditions are satisfied on the plane between the z mesh positions k m and k m + 1 
or z = (zk m + z k m +l )/2. For convenience, we choose this plane as 

Z = Zlr ■ . .„ = 0 
K m+l/2 
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The boundary conditions in equation (28) modify the V\ zz derivative in the difference 
equation. At the point i j,k m , we have 


2 


^ijk m - ^ijk m -1 _ 

Zk m +1 ' z k m -l 

^lz 

Zk m ' Zk m 4 . 


z z k m +l/2 


( 29 ) 


Applying the boundary conditions, equation (28), leads to 

^lzz = 2h l b zk m F i/ L) ' 2a zk m ^ijk m ‘ ^ijk m -l^ (30) 

where 

F ..(L) = F( L )( Xi , yj ) = f j x (L) (xi,yj) + iwf ! (L) (Xi, yj ) 

and 



+ 1 ' z k 


m 


Similarly, for the value of V\ zz at z = zk m +l> we have 

2 pjk m +2-njk m + l I-! 

zz ~ ( z k m +2 ‘ 2 k m )l_ Z k m + 2 ‘ z k m +! ^'"1 J 

z = z k m +l/2 

m (31) 

^lzz = 2b zk m +l (^iik m +2 ' ^ijk m + 1 ) * 2b l a zk m + 1 F i/ ^ 

The. superscripts L and U denote upper and lower surface values, respectively. 

For k = k m , the difference equation corresponding to equation (24) is 

a zk m ^ijk m -l * ( a yj + b yj + a zk m + E 1 + E 2 " q ijk m / 2 Mjk m 
= - E 2^-ljk m - a yj <# ’ij-lk m ' byj^ij+lk m - h, < 32 > 
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where 

E l =C i U i+‘/2jk m - iWC li/ e 
E 2 = d i u i- 1 /2jk m + ia,d li / £ 

Similarly, for k = k m + 1 , we have 

-(a y j + byj + + Ej + E 2 - Qijk m +l/^^ijk m +l + b zk m +l^ijk m +2 

(33) 

= - E l^i+ljk m +l ‘ E 2^i-ljk m +l * a yj^ij-lk m +l ’ b yj^ij+lk m +l + h l a zk m +l F ij (U) 

where 

E 1 = C i u i+V 2 jk m +1 ' iwc li/ e 
E 2 = d i u i->/2dk m +l + 

Analogous equations are easily written down for hyperbolic points by comparing equations (24) 
and (27) and revising (32) and (33) accordingly. 

Boundary Conditions on the Wake 

The pressure and normal velocity component must be continuous across the vortex sheet shed 
from the wing trailing edge. At z = (zk m + zk m +l)/2, these conditions become 

*>f x + iwcpj= <p'i x + ioj^'i (34) 

for pressure continuity, and 

Az = z 
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for normal flow continuity. The superscripts + and - designate limits as the sheet is approached from 
above and below, respectively. Equation (34) may be written 

(d/dx)(Aipj) + icjAipj = 0 (35) 

where Ai^j = *P\ Along each line y = constant representing the linearized streamline from the 

wing trailing edge, equation (35) may be integrated to give 


A<pjj = Aipj i + ! j e~ ia> ( x i" x i j + 1 ^ (36) 

where A^ij+ij is the jump in the potential at the station yj for the first mesh point x location 
beyond the trailing edge on the plane z = 0. This value is determined to satisfy the Kutta condition, 
which requires that the pressure on both sides of the trailing edge be equal. Let ij be the index for 
the value of x on the trailing edge. Then the condition that the pressure jump vanish at the trailing 
edge at the station yj becomes 

c li 1 ( A ^i 1 +ij- A ^i 1 j) + d li 1 ( A ^i 1 j' A ^i 1 -lj ) + iwA ^i 1 j = 0 (37) 


Solving for A<pj^ + jj yields ' - 

A ^i 1 + lj = A ^i 1 j [1 “( d lij +i ^)/c li ] + (d lil /c lil )A^ il . lj (38) 

To satisfy continuity of ipj z across the vortex sheet, although may be discontinuous, we 
write for ^j 1z at z = zk m + 1/2 = 0, 


.^z . (<p ijk m +l *. A ^j -^ijk m .)/( z k m +l * z k m ) 


Hence, for k = k m + 1, the derivative <£\ zz becomes 


*lzz " 2a zk m +l ^ijk m ' 2 ( a zk m +l + b zk m +l^ijk m +l 
+ 2b zk m +l ^ijk m +2 + 2a zk m +l A *ij 


(39) 
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Similarly, for k = k m , we have 


^lzz 2a zk m ^ijk m -l - 2(a zk m + b zk m Mjk m + 2b zk m ^ijk m +l ' 2b zk m ^ij 


(40) 


By comparing equations (39) and (40) with equation (22) we see that, in the difference 
equations (24) and (27) for k = k m , the right-hand side has the additional term 


b z k m Aip ij 


(41) 


while in the difference equation for k = k m + 1 , the right-hand side has the additional term 


■ a zkm+ lAl/, ij 


(42) 


By employing this form of the difference equations, we implicitly apply the condition of continuity 
of i/3j z at the plane z = 0. 


Trnasformation of the Differential Equation for Swept Wings 

Since a finer mesh is desirable in the vicinity of the leading edge, it is convenient to transform 
the coordinates so that the swept wing is rectangular with the leading edge as one of the coordinate 
axes. Let the leading edge be defined by x = y sin a or y = x/sin a. Then we introduce coordinates 

% = x - y sin a 

v = v (43) 


? = z 


The derivatives with respect to x and y in terms of the new variables then become 

^lx = ^U 

^ly “ ^lrj ' sbl a ^1| =i ^1tj'^1£ 
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and 


*izz _ *ltt 

^lyy = *lr)V ' 2X *1$t? + 


where 


X = sin a 


Substitution into equation ( J 7) then yields 


[vy>j£ - 2X^1^ - (2ico/e)v?j ] | + + qv?j - 0 (44) 


where 


v = u + sin^ a = u + X 2 


We consider the second-order operator of equation (44) and apply Hormander’s theory (ref. 1 6) 
to find under what conditions it is hyperbolic and to define the cone of influence of a point distur- 
bance in %,T},t Space, The operator is 

P(D) = vDj 2 - 2 XDjD 2 + D 2 2 + D 3 2 

where D n = i9/9x n with x i = £, *-2 = ^ x 3 = ?• According to Hormander’s definition (ref. 16), the 
operator is hyperbolic in the £ direction if 

P(1,0, 0) = v = u + X 2 =£0 

and if 

P (D + rej) 


has only negative roots for D real. Here e^ is a unit vector in the £ direction. The roots of 
P(D + rej) are found from 


v(Dj + r) 2 - 2X|D.j + r)D 2 + D 2 2 + D-^ 2 - 0 
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Solving for r yields 


r = -Dj + (XD 2 ± ^-uD 2 2 -vD 3 2 )/v < 0 

Since r must be real for D real we see that u and v must be negative or u < -X 4 . The condition 
r < 0 becomes 

D 1 2 v-2XD 1 D 2 + D 2 2 + D 3 2 >0 (45) 

The cone of support of the solution in physical space must also satisfy 

£Dj + > 0 

for Dj , D2, and D3 satisfying equation (45). Both regions are convex. The boundary of the region 
is defined by the equality of the preceding two equations. Eliminating D3 from equation (45) by 
means of the preceding equation yields 

F(/Jj ) - ' 2 a 2 v + £ 2 ) - 2& x (Xf 2 - + (r? 2 + ? 2 ) = 0 

with jS j = Dj/D 2 . 

We now find the envelope in the physical space by finding the value of j3j for which 

F(^ 1 ) = F'(|3 1 ) = 0 

The equation for the envelope is found to be 

f 2 [u(? 2 + r? 2 ) + (i?X - £) 2 ] =0 

Let £ 2 + £ 2 = r 2 and 17 = r cos 0, then we have, after solving for r, 

r/£ = (X cos 0 ± vpu)/(u + X 2 cos 2 0) 
or 

r-J-u/l; = (i^ cos 0 ± l)/(^j 2 cos 2 0 - 1) 

where 

l ; 1 =X/ s / : u<l 
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Taking the minus sign gives us 


«%/*>/£ = 1/(1 - v \ cos 0) 

This describes the cone of influence of the point 0,0,0 downstream of the point. Taking the plus 
sign yields 

r s/-u/£ = 1/(1 + cos d) 


This equation represents the boundary of the upstream region of influence on the point = 0,0,0. 

For the rectangular wing, X (or v j ) goes to zero and the equations reduce to the usual circular 
Mach cones for the regions of influence. For -u < X, the cone cross sections for f = constant are 
ellipses, while for the limiting case = 1 , or X = -u, they are parabolas and, of course, are not 
closed. The cross sections in the t? = 0 and f = 0 planes are bounded by two intersecting straight lines. 

Since for u < -X^ or <pq x > (K + sin^a )/(7 + 1 ), the transformed differential equation is hyper- 
bolic, we employ backward differences in the £ variable. Central differences are used when 


^0 x < (K + sin 2 a)/(7 + 1) 


(46) 


Difference Form of the Differential Equation for the Swept Wing 

We first establish a second-order difference for . By comparison with the second-order dif- 
ference for ip] X) we write 

*ln = c lj^ij+lk “ ^ijk> + d ij^ijk ‘ (47) 

v 

where cy and djj are given by the equations following equation (19) with r? replacing x and j 
replacing i. 

Applying a second-order central difference operator for | to equation (47) yields 

^1£tj = c lj 0lM+lj+lk ‘ %Hk> + d lMj+lk ' ^i-lj+lk^ 

+ ( d lj - c lj )[c li («p i+ ijk - ^ij k > + d iMjk - ^i-ljk)] (48 ^ 

- d 1 j [c i i(^ i+ lj-lk - ^ij-lk> + d l Mj- 1 k ' ^i- 1 j- 1 k^ 
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Then 


-2**!*, - 2(ico/eV 1 | = ^XcjjlcjiCvJj+jj+jk - ^ ij+lk ) 


+ d li^ij+lk-^i-ij+ik) ] 


- 2[ic o/e + X(djj - cjj)] [ c iiOPi+ijk " ^jjk) + d li^ijk ' ^i-ljk^ 


+ 2Xd lj [c li (v? i+1 j. lk - ^ij-ik) + d li^ij-lk ' ^i-lj-lk^ 


For the other derivatives we have 

lw>]|] | - 2cjVj + ^ jk^i+ljk " ^ijk) * 2d i v i- L /z jk^ijk ' ^i-ljk) 

^1777? = 2 b 7 y^ij+lk ‘ ^ijk) " ^ a 7 ?j^ijk ' ^ij-lk) 


( 49 ) 


(50) 

(51) 


*1 ff “ 2b ?k ( ^ijk+l ‘ ^ijk^ ’ 2a fk^ijk ' ^ijk-1 ) ( 52 ) 

where a^j, b^j, aj- k , and b^ k are defined by equation (23) with appropriate variables and subscripts. 
Substituting equations (49) through (52) into differential equation (44) and collecting terms yields 

a £k^ijk-l " + b T?j + a ?k + b ?k + E 2 + E 3 " 4ijk/ 2 Mjk 

+ b fk^ijk+l = ' E 4^ij+lk ‘ E 5^ij-lk ' E 2^i+ljk 

( 53 ) 

' E 3^i- ljk + Xc 1 jl c li^i+lj+lk * d l i^i- 1 j+ 1 k 1 

- Xdjjtcjj^j+jj.jk - djj^.jj.jk] 

where 

Ej = icu/e + X(dy - cy) 

E 2 = c i v i+ jk - E l c li 

E 3 = d i v i- 1/2 jk + E l d li (54) 

E 4 = ^ c lj( c li " d lP + b 7jj 

E 5 = - Xd 1 j(c u - d n> + 
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Equation (53) is the difference equation in swept wing coordinates equivalent to equation (24) for 
elliptic points. 

For hyperbolic points, we have v < 0 or u < - X . For the £ derivatives we qse backward 
differencing. Thus 


[v*^] £ = 2c i _ 1 v i _ 1/2jk (v> ijk - <Pi_i jk ) (55) 

' 2d i-l v i-3/2jk^i-ljk ' ^i-2jk> 

-2X^1 ^ - 2i6J(pj£/e= -2Acjt c 2i(^jj + ik - ^i-lj+lk^ ' d li-l^i-lj+lk 

" ^i-2j+lk^ * 2E 1 t c 2i^ijk ' ^i-ljk^ ' d li-l^i-ljk ' ^i-2jk^ 

+ 2Xdj [ c 2i (vP i j. lk " ^i-lj-lk) " d li-l ^i- 1 j- 1 k ' ^i-2j- 1 k^ - 


where c 2 j is defined by equation (26) with £ replacing x. Substituting equations (5 1), (52), (55), 
and (56) into the differential equation (44) yields 

Wijk-1 ' ( a r?j + b 7?j + a fk + b £k + e 6 + %‘k/ 2 ) ^ijk 
+ *Wijk+l = E 8^ij+lk ' E 9^ij-1 k + < E 7 ‘ E 6M-ljk 

(57) 


' E 10^i-lj+lk + E 1 l^i-lj-lk + E 7^i-2jk 
+ Xc 1 j d 1 i- 1 ^i-2j+ 1 k ' Xd 1 j d 1 i- 1 ^i-2j- 1 k 


where 

E 6 = c 2i E l ‘ c i-l v i-l/2jk 
E 7 = d li-1 E 1 - d i-l v i-3/2jk 


E 8 _ Xc lj c 2i ■ \} 


e 9 - Xd ij c 2 i + a rjj 

E 10 = Xc lj( c 2i + d li-l) 


( 58 ) 


E 1 1 _Xd l/ c 2i + d li-l ) 
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and E j is defined in equation (54). This is the difference equation for hyperbolic points in swept 
wing coordinates equivalent to equation (27). 

Root Chord Boundary Conditions for the Swept Wing 

To limit the field to one-half the span for rj > 0, we require the symmetry condition 

ip ly = 0 at y = 0 

or in swept wing coordinates 

^Itj + X< ^1£ = 0 

Since this is true along tj = 0 for all f, then we also have 

* 1t?£ + X ^U£ = 0 

This relation conveniently enables us to eliminate the cross-derivative term Thus, on the line 
of symmetry, the terms in the differential equation (44) involving the f derivatives becomes 

[vipj g - 2Xipj ^ - (2ico/eV j ] = [(v + 2X 2 )<£j £ - (2icu/e)i/?j ] (59) 

Let r? = 0 lie on the j = 1 horizontal mesh line. Then for points j = 0, 1, and 2, the difference form 
of the boundary condition is 

*1 r? + Xv? l£ = ( ^2k ' ^i0k^2 ' X[c lM+l lk '^ilk^ , 

+ d lMlk'^i-llk)l =0 

We choose for convenience tjq = -77 j . Then for elliptic points, we solve for and obtain 

^iOk = 9 ? i2k - X (ri2 + *?l)lciM+l lk ’ ^ilk> + d lMlk ' *i-l lk>l 
Similarly, for hyperbolic points 

^iOk = ^i2k ‘ X ^2 + *?lH c 2Mlk * ^i-1 lk)' d li-l( ^i-1 lk ' ^i-21k>l 
Also for j = 1 , we have 

T7T? = 2 ^i2k ‘ ^i 1 k)/^2 2 ' ^ 1 2 ) ' ^i 1 k * ^iOk^ 1 0?2 + ^ 1 ) 
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Eliminating i/>jQ k for elliptic, points yields 

* 1 T?T? = ^i2k - *i 1 k)/ 7 ? 1 (*2 - ** 1 ) - W 7 ? 1 ) [ c 1 M+ llk'^ilk) 

+ d lMlk'^i-llk)l 
Similarly, for hyperbolic points 

* 1 r?T? = 2k - Vi 1 kV 7 ? 1 (*?2 • ^ 1 > - W 7 ? 1 > [ ' c 2i^i 1 k ■ *i - 1 1 k> 

- d 1 i- 1 (^i- Ilk" ^i-2 1 1 

For j = 1 the £ derivative for elliptic points is 

[(v + 2\2 )i£j£ - (2ico/e)<pj ] £ = 2E| 2 (y>j+i - ^j]k) 

-2Ei 3 (^iik-sPi-iik) 

where 

E 12 = c i^ v i+’/2lk + 2X ^) - icocjj/e 
E 1 3 = d i( v i-%lk + 2x2 ) + icodjj/e 

and for hyperbolic points 

[(v + 2A 2 ) v ? 1 £ - aioj/eVj 1 £ = 2E 14 (<p ilk - ip M lk ) 

- 2E 1 5^i-l lk - ^i-21k) 

where 

Ej4 = Ci-i (Vj.x^ i k + 2X 2 ) - iojc 2 i/e 
E 15 = d i-l( y i-3/2 lk + 2x2 ) - iwd li-l/ e 


( 60 ) 


(61) 


(62) 


Substituting equations (62), (60), and (52) into the differential equation yields the following differ- 
ence equation for j = 1 . 

a ?k^ilk-l ' [b r?l +a fk + b fk + E 12 + E 13 + X ( d li- c lP/ 2? 7l 

^i 1 k/ 2 ] ^i 1 k + b r k^i 1 k+ 1 = ‘ b r? 1 ^i2k 
+ (Xc li /2r? 1 -E 12 V i+llk -(\d li /2r? 1 + E 13 V i . nk 



where 


b T?! = 1/2 i?i(t7 2 - 77J) 


Similarly, for hyperbolic points, we have 

a fk^ilk-l ' ( b 77j + a ?k + b fk + Xc 2i/ 2l? l ' E 14Mlk 

+ b rk^ilk+l = ' b r? 1 ^i2k + l E 14 + E 15 '( x / 2T ?l>( c 2i 

+ d li-l)]^i-llk' (E 15 ' M li-l/ 2T ?lM-21k- 

These last two equations replace equations (53), (54), (57), and (58) for j = 1 points in the mesh, 
which are the points on the plane of symmetry for the flow. 


Wing and Wake Boundary Conditions for the Swept Wing 


The boundary conditions on the wing and on the wake involve only the £f term, which is 
unaffected by the coordinate transformation. The difference equations for k = k m and k m + 1 are 
found by modifying some of the terms in equations (53) and (57). For k = k jn , b£k m deleted 
from the coefficient of ^ijk m ) tbe ^ijk m +l term is omitted, and 


-h 



..(L) 

ij 


is added to the right-hand side. Similarly, for k = k m + 1 , the <£ijk m term is deleted along with 
a fk +1 from the coefficient of + j and 

l ’i a ?k m +i F ij (U) 


is added to the right-hand side. 

The wake boundary conditions for rectangular wings also may be carried over into the swept 
wing coordinate system with little change, since 

*lx = *l£ 


The second equation following equation (35), and equation (38), can be applied with only a change 
in notation. For k = k m , we add to the right-hand side of equations (53) and (57) the additional term 
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b rk m 

and for k = k m + 1 , the additional term - (63) 

- a fk m +l A ^ij 


Far-Field Boundary Conditions— Three-Dimensional Flow 


In order to solve the difference equations in the flow field, we must prescribe boundary condi- 
tions on the exterior boundary of the mesh. This is usually found in terms of the velocity potential. 
This boundary information is not known a priori but must be computed and periodically corrected 
in the iteration process used to find the solution. For this purpose we follow an analysis similar to 
that of Klunker (ref. 8) for the steady-state equation and write the differential equation in the form 


LCspj ) = Ki/>j xx - (2iw/eV lx + ip lyy + f\ zz + (oj 2 /eVj 
= (T + n(v?0x^lx)x + - 1 Vl^Oxx 


(64) 


To convert the solution to the swept wing coordinates, we merely substitute x = £ + X 77 and y = rj 
into the solution of equation (64). For our purposes, it is convenient to introduce new variables in 
terms of x and y. Let xj = x, yj = y V^”and Zj = zVlCThen 

f2 

^lxj + ^lyjYj + ^IzjZj 


L ^l) = ^i Xl x 


1 X 1 


2icoM‘ 


1 -M" 


co 2 M 2 


1 -M“ 


1*1 


_ 7 + 


i co(7 - 1 ) 


K - (*0x]*lx]'xj + K *l*0x 


1 X 1 


We derive an integral equation for ipj in terms of boundary values by applying Green’s theorem in 
the form 


f. 


WLOpp-vJjUMdv = 




^ln 


2icoM" 


1 - 


<pj^n x Ids' 


(65) 


where the volume v is bounded by the surfaces S = Sp + + S w + S v (see fig. 18). Here Sp is the 

sphere about the specific point p. The surface is the surface of discontinuity representing shocks. 
The quantities S w and S y represent the wing and trailing vortex sheet, respectively, as shown in 
figure 1 8. The quantity n is the normal directed into the fluid from the surfaces. 

The fundamental solution \p is a function of the vector Xp - x\ where Xp is the vector refer- 
ence point and x' is the point vector of integration. 
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To find the fundamental solution \l/, we let 


icoM 2 xi/(l -M 2 ) 

^ “ e 


With r - Vx] 2 + yj 2 + z j 2 , the differential equation, L(i p) = 0, becomes 

^ o " + 2\}/q/t + \ 1 2 ^q = 0 


where Xj 
leads to 


= coM/(l - M^). The solution of this equation. yielding outgoing waves is e 


-iX, 


r 

/r and 


iXi(Mxi-r) 

\ // = e /r 


( 66 ) 


for the fundamental solution. We let the radius of the sphere about the reference point Xp = (Xj , 
Yj,Zj) shrink to 0. Since \p n = i//^, then the contribution of the surface integral on the right-hand 
side of equation (65) about Sp yields 


(67) 
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- 47 r^xj, yj, zj) 



Since L(i//) = 0, the volume integral in Green’s theorem becomes 



[i//L(y3 1 )-(£ 1 L(i//)]dv' 


) 


= ' K l 


[(7 + UMixjVOxjOxj' + M?- 


where the primes on the variables denote the running variables of integration. Integrating the first 
term by parts with respect to x j ' and noting that a jump occurs on the shock surface yields 


(7+1) 
K 


i 



^Ox^lxj'l costs' 


where cos 0 j is the inchnation of the normal to the shock surface with respect to the xj axis. The 
square brackets denote the jump in the quantity across the surface. For the integral over the shock 
surface, we obtain 



where cos an <3 cos @3 are the yj and zj direction cosines of the normal to the shock surface, 
respectively. 


For the thin planar wing, m then the surface integral over the wing reduces to 



/•y t fX&O 

(M n -^ln) ds '“/ / , ( A ^l^z/ 

•'-yt •'xgCyj ) 1 


^ A V’l Zl ».)dxi'dy 1 < 


(69) 


where Xg and x t designate the X] coordinates for leading and trailing edges, respectively; y t is the 
value of yj , for the tip of the wing, and A denotes the difference between upper and lower values 
of the quantity on the wing, which is approximated by its projection in the plane zj = 0. 

Similarly, the integral over the trailing vortex sheet, since v?j n is continuous, yields 



(y’l'/'n' ^i n )ds' 



A< ^l^Zj dxp dy^' 


(70) 
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Now, continuity of pressure across the vortex sheet requires that 


(Av>i) X] ' + icoAip] = 0 


Solving the differential equation yields 


A«pj(xi ,yj ) = A^ t (yj )e 1 1 1 


( 71 ) 


(72) 


where A^ is the jump in i pj at the trailing edge of the wing. Substituting equation (72) into equation 
(70) yields for the wake integral 



icox 


{y\ ) 


A<^ t (y i )dyj 


r 

JxJv 1 


-1CJX 


x t(y i ) 


V z, dx l 


(73) 


We have now found the integrals over the volume and over all surfaces of figure 1 8. Substituting 
equations (67) through (73) into equation (65) and combining integrals over the discontinuous sur- 
face S^ yields 





2icoM 2 
1 -M 2 


[t^l 1 n x ) COS0] 


+ cos(9 2 + cos ^ 3 1 ds' 


(74) 


This can be shown to vanish identically for shocks by using the divergence theorm. The partial differ- 
ential equation can be written 


F 1 x j +F 2y 1 +F 3Z] + Q^1 -° 


(75) 


where 



I+i 

K 


^Oxj^l x j 


2icoM; 


1 -M 


2^1 


F 2"^l yi F 3 = ^l Zl 

^ ioo 2 M 2 . , 

• q= ~Tm 2 ■ ic °(y ■ )^ 0x j x j • 
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Consider the integral over a small rectangle about the shock as shown in figure 1 9. Applying the 
integral of the differential equation (75) over this small volume and letting the thickness At-K) 
yields from the divergence theorem 

(V * F + Q)dv f =*A s[F • n] = 0. (76) 


hm' f 

At-*-0 *' v 


The quantity [F • n] is seen to be identical to the quantity enclosed by j } in equation ( 74 ), since 
n «= (cos 6 j , cos 6 2 , cos 6 3 ). This demonstrates that the integral in equation (74) is identically 0. Thus, 
the Green’s theorem of equation (65) becomes 


/*yt/* x t(yi') 

= £ , z'ldx^dyj' 

47 r 4 t j h(y\ ) 1 


*1 


J -y t ‘ /x t^ y l ) 


+ &rkf y ^ y + D^i 



(77) 


Figure 1 9.— Shoe k surface and volume used for divergence theorem, integral 
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The wake integral can be integrated with respect to x j ' for the combination <pj- x + icoipj , which 
may be called the pressure potential. We consider 


/oo 

-icoxi' , , , , 

e +iw^ Zl '] dx l 


since 9/9xj = -9/9xj\ we have 


w 


_ f°° 9 . -iwx/ , . -iwx t , , 

-J aTj' ( e 1 ^ Zl ') dx i = ~ e * Zl < x r x i »yi-yi > z i>- 


(78) 


Since the wake integral is very slowly convergent, it is more convenient to express the upstream and 
downstream boundary conditions in terms of <£] x + icoi^j rather than <pj. The pressure function from 
equation (77) then becomes 


, nt r x t(yi'y 

<Pl x + koy>i =4- / / [A^> ]X 2 ’-xAviz 'Idxj'dyj' 

-yt *Vyi ) 


1 /* yt icox t (yi') A , . 

-ml e A ^ t (yi )^ x t(yi )- x i>yi -yi > z i) d yi' (79) 


-y t 


4^k/*( ( 7 + ^Oxj' Vlxj'Xxj' - X<P 0x ^ Xi - ) dv' 


where 


X= +ico^ 


Formulation of the Problem of the Harmonically Oscillating Airfoil 

To test the method developed in the preceding sections, we formulate the relations for the har- 
monically oscillating airfoil, since they are simpler and will require less computing time. The equations 
for the oscillating wing represent the wing as lying in the plane z = 0 with the direction of the undis- 
turbed flow in the positive x direction. Relations were derived for the method of relaxation of the 
difference equation by solving for the values of the potential along lines normal to the z = 0 plane 
corresponding to fixed values i and j points in the mesh. To follow the conventional notation for 
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airfoil problems, we now consider the airfoil cross section as lying in the x,y plane. Since there is no 
variation with z then the small perturbation differential equation (17) for the airfoil becomes 

[u<p lx -(2ko/eVj] x + <£j yy + q<P] =0 (80) 


where 


u = K - (7 + l)<^0x 
q = co 2 /e - iw(7 - l)<po X x 


(81) 


By representing the airfoil in the x,y plane we effectively interchanged y and z in the derivation of 
the fundamental equations. The unsteady boundary conditions on the airfoil corresponding to equa- 
tion (28) then become 

ip ly = f 1 '(x) + icof 1 (x) (82) 


on y = 0, where yQ = 5 fg(x) + 8 f j(x)e lwt describes the oscillating airfoil shape. The boundary con- 
ditions of equation (34) become continuous across y = 0 behind the airfoil, and equations (71) 
and (72) for continuity of pressure lead to . 


a a -icoCx; - x t ) 

Aipj = Ai£ t e l Xj>x t , 


(83) 


where is the jump in the potential at the first mesh point downstream of the trailing edge. Since 
we have dropped variation with respect to z the difference form of the differential equation (80) 
is found from equation (24) with a z ^ = b z j c = 0. For interior points in the flow field we then have 


a j^ij-l ‘ ( a j + bj + Ej + E 2 - q ij /2)p i j + bj <^ j+1 = - E^ i+1 j - E 2 </> Mj , 

where 

E 1 = c i u i+l/2j - E 2 = +icod u /e. 


(84) 


and we have dropped the subscript y on the coefficients a and b. 


The boundary conditions on the airfoil are satisfied at y = 0, the midpoint between y = yj and 
yj + j . From equations (32) and (33) the difference equations corresponding to equation (84) ^r 


Jm 

J = j m 
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and j m + 1 become 


a j ^ii -1 " ( a j + Ei + E 2 - q- /2W-: = ' E l^i+li ‘ E 2^-li " h 1 b i f/ E \ 

J m JJ m J m 1 z ] Jm m 1 IJ m zn IJ m 1 J m 1 


( 85 ) 



and 


where 


Here 



+1 +Eo 
m J 


+ E4 - 


q ij m +l 


m 

2 . 


Hm + 1 b jm +1 ^m +2 


-- E 3^i+lj m +l - E 4 ^i_i j m +i +h 1 a jm +l F i (U) , 


( 86 ) 


E 1 = c i u i+l/2j m - i0JC U / e > E 2 = d i u i-l/2j m + icod li/ e 

E 3 = c i u i+l/2j m +l - iwc li / e > and E 4 = d i u i-l/2j m +l + iwd li/ e 


(87) 


F 4 = f j '(Xj) + icof 1 (X)). 


Downstream of the wake the boundary condition of continuity of pressure and normal velocity 
component require that, for j = j m , equation (84) holds with the additional term on the right-hand side 

b; Aipj 
J m 


and for j = j m + 1 , the right-hand side of equation (86) has the additional term 



Some of the coefficients of the difference equations for the harmonic motion of an airfoil in 
transonic flow require the steady-state solution. For this we use a program called TSONIC or TEA- 
330, which was developed by J. A. Krupp (ref. 4) and computes the solution of equation (14) for 
the transonic flow over lifting airfoils. The output of this program is an array of values of the pertur- 
bation potential y>g corresponding to the mesh points of the rectangular grid about the airfoil. 

Since the coefficients in the differential equation depend only on <^q x and ^Oxx 5 ^ conven i en t at 
the outset to establish values of u = K - (7 + 1 )i/>q x at each mesh point, since <^q x occurs only in this 
quantity in the equation. For convenience of coding on the computer we write 

u i-l/2j "*• «ij = K - (7 + D(^0ij -^Oi-lj>/< x i - x i-l) 

and (88) 

Uj+l/2j u i+lj = K - (7 + l)(<Poj + ij - ¥? 0 ij)/( x i+i - Xj). 
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Only the quantity q contains the term </>Qxx : 


q 


= co 2 /e ■ 


ico(7- 


^Oxx- 


From equation (88) we have 


(^Oij ‘ ^Oi-lj) = ( x i - x i- 1 X K - U ij)/(T + 1 ) 


Using this relation in the difference form for <pq xx 

^Oxx = c i%i+ 1 j ' Y’Oij) * d i^0ij ’ ^Oi-lp 

= [(x i+1 -x i )c i (K-u i+1 j)-(x i -x i . 1 )d i (K-u i j)]/( 7 + 1). 


( 89 ) 


(90) 


where cj and dj are defined following equation (19). Simplifying equation (90) leads to 

*>oxx = ( u ij ' u i-lj)/(7 + IKxi+i - x i_i). 


Finally, for q y we write 

q ij = co 2 /e+ic 3i (u i+lj -u ij ) (91) 

where 

c 3i = w(7 - 1 )/( 7 + 1 )(x i+ j - Xj. j ). (92) 

For the volume integrals in the far-field integral relation, we require ipg x at the points ij. For this we 
write 


*0x " c lM)i+lj ' *0iP + d li ( ^0ij ' ^Oi-lj) 


which, with the aid of equation (89), simplifies to 


(7+ 1 Vox = K ’ c 4i u i+lj " d 4i u ij 

where 


(93) 
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c 4i = <- x i - x i-l)/Cx i+ i - x i_i> 


( 94 ) 


d 4i = ( x i+l ‘ x i)/( x i+l ' x i-l>- 


Equations (91) through (94) provide all the relations for the steady-state solution needed in the dif- 
ference formulation of the harmonically oscillating airfoil. We have replaced the v?Qjj array by the 
uy array defined by equation ( 88 ). 


Evaluation of the Jump in on the Airfoil 


Since the values of <p j are computed only at mesh points, we must interpolate these values to 

the limit y = 0 in order to compute the values of Aipj required in the wake boundary conditions. 

The limiting value of from above the airfoil and from below are found by Taylor’s 

expansion about the points y = y: and y = y; + i. We obtain 

J m J m 1 


<Pl (L) = V>ij +(h/2V ly l. +(h 2 /8)^ lyy |. +••• 
m Jm J m 

^! (U) = ^ij m +l - ( h /2V ly | +(h 2 /8)^ lyy , +••• 

111 J m 1 J m 1 


where h = Yj m +1 ■ Yj m « We use a second-order difference form for <pjy. Following the formula 
developed for <p\ x , we write for 


% Jm +1 =C3(< ’ ijm+2 ^ ijm+l) + d3<, ’ iim+1 ~ Vl<U)> (97) 

*lylj = °4<*’1 <L) • + d 4 < f , ij m (98) 


By comparison with cy and dy following equation (20), we obtain 



C 3 = l/hsj( 2 sj + 1 ) 

d 3 = 4s 1 /h(2s 1 + 1) 


c 4 = 4s2/h(2s2 + 1 ) 

d^ = l/hs 2 ( 2 s 2 + 1 ) 

where 

s ' =(y im«- y i m + i )/h 

and s 2 = (y; -y; . 

^ J m J m 


(99) 
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The second derivatives at y = and y = y: , i become 

J m Jm 1 


J m 


-1 


” lyy ij m = yj + i- y j -i 


t F i 


.(L). 


m 


m m 


y i ’ y i 
j m J 


(100) 


m 


= 2Fj( L )/h(s 2 + 1) - 2(^j -^ ij . 1 )/h 2 s 2 (s 2 + 1) 

m m 

*lyyl. , 2 77, (y ij +2 - +l)-2F i <U>/h (S i + 1) 

J m 1 h s l< s l + n m m 


(101) 


Substituting equations (97) through (101) into equations (95) and (96) and solving for and 

i leads to the following relations after some simplification: 


^ 1 ( 0 = + (^y - ^|j _j)/4s 2 (s 2 + 1) + h(2s 2 + 1)F^ L )/4(s 2 + 1) 

m m m 


( 102 ) 


and 


^](U) = (/? „ +J -(v?y + 2 - <^ij +l)/4s 1 (s 1 + l)-h(2s! + l)F i ^ U ^/4(s 1 + 1) 
m m m. 

Finally, the jump in <p . j is given by 

Atp = n (U) . ^j(L) = ^ +1 - <py - C sl (^j +2 - +1 ) 

m 1 1 mm 

- c s2<^j m - .i)-(d sl F i (U > + d s2 F i (L, ) 

111 m 

where 


(103) 


(104) 


c sl = l/ 4s l( s l + 1) c s2 = 1 /4s 2 ( s 2 + 1) 

(105 

d sl = h(2sj + l)/4(sj + 1) d s2 = h(2s 2 + l)/4(s 2 + 1) 

These formulas hold only for values of i corresponding to x values on the airfoil and are accu- 
rate to the second order in mesh size. The jump in <pj is required for updating the far-field boundary 
conditions as well as satisfying the Kutta condition on the airfoil trailing edge. 
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Far-Field Boundary Conditions— Airfoil 


A procedure similar to that used for three-dimensional wings is employed to obtain an integral 
relation for the far-field boundary condition for the two-dimensional airfoil. We write for in 
place of equation (64) 


2icoM 


w 2 M 2 


+ ^iv,yi + T^5 Vl 


7 + 1 , v , ico (7 -. 1 ) 

K ^Oxj^lxj^X] + K ^I^OxjXj 


(106) 


For the fundamental solution t// of the differential equation L(<//) = 0, the function i//g(r) takes 
the form 

Pq” + 'Po'/r + 2 ^o = 0 


where r = Vxj 2 + yj 2 and Xj = coM/(l - M 2 ). The solution to this equation yielding outgoing 
waves with the source at the origin is HQ^ 2 -*(Xjr). The fundamental solution then is 

iXiMxi 

* = e 1 H 0 ^ 2 ^(Xjr) 

where Hq^ 2 \z) is the Hankel function of the second type. For small r 


, 1 \ iXi Mr cos 6 . ... 

P “ -(2i/rr)e 1 log(Xjr/2). 


Since \p n = p r the integration of the surface integral of equation (65) over the circle about the 
reference point xj,yj becomes, on letting the radius of the circle go to zero, 


- 4 iip^xj, yj). 


(108) 


The evaluation of the remaining integrals for the two-dimensional flow follows with little change 
from the three-dimensional analysis. Noting that the factor multiplying <p j is -4i insteady of -47T, 
we obtain easily the following equations for the airfoil equivalent to equations (77): 


^l - (l/4i)^* (Aip j \py j' - i//Ai/>iyj’)dxi + (A<Pf./4i) ^ e ^ ^ ^y^'dxj' 
+ (l/4iKj^ [(7 + 1 )^xj'^0xj' ^lxj' - ico(7 - 1 VQxj'x^I tMdxj'dyj' 


(109) 
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The quantities in the first two integrals are evaluated at yj =0. Since 9/9xj' = -9/9xj, 9/9y j ' = 
-9/9yj', equation (109) takes the more convenient form 


<Pl( x l,yi) = G/4)^ (A^ j ^ y j + ’/ /A ^iy 1 ')dx 1 '' +(iA<p t /4j ^ 

+ (i/4K )J [(y+ 1 )'/' Xl 'flo Xl Vix 1 ' + iw (T- 1 )‘fl3x 1 'x 1 , ^l ^dxj'dyj' 


^ yi dx l' 


( 110 ) 


As shown for the planar wing, the integration of the second integral can be performed for the com- 
bination +ico<^j. Thus 


^1X] +iw ^i =mf (Asp iXyj +xAv? 1 y i ')dx 1 ' + (\A<p t /4)\jj y ^ 
W [(<y + ’^xj^Oxj'^Ixj' +i ^(T- IVoxj'xj'^lXldxj'dyj' 


dll) 


The quantity X = ’/'xj + Iuji// I s a function of the arguments xj - Xj \ and yj in the first integral and 
of Xi - x/ and y j - y j ' in the volume integral. In terms of x and y, we obtain 



+ xA<p ly ')dx' + (iA<p t /4 s/K)t// y 


(i/4 s/K ) j [(7 + DXx'POx'^lx' +ioj(7 ‘ 1 )'P0x'x'^l x]dx ' dy ' 

Jv 


( 112 ) 


where 


* = e iMX l X H 0 < 2 >(X, 

i>y - -AjKy e lM ^l >; H| ( - ‘ (X ^ -''x - + Ky-J/ ''X^ + Ky~ 
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x = ic O^K 1 - M 2 ) - Xj X e lXl Mx Hj < 2 )(\ j 7x 2 + Ky 2 ) / J x 2 + Ky 2 
X y = icov// y /(l - M 2 ) + XjKxy e^ 1 ^ [2 Hj( 2 >(Xj Vx 2 + Ky 2 ) / Vx 2 + Ky 2 

- H 0 ( 2 )(X j 7x] 2 + Ky 2 )] /(x 2 + Ky 2 ). 


Equation (112) holds everywhere in the flow field. Since we are only interested in the far-field 
boundary values, we approximate x - x' in x and Xy of the integral over the wing by x. The terms x 
and Xy can be taken outside the integral, and the first term in equation (112) then takes the form 


(i/4 7 k) tx y ^ j dx- + xjf 


A<p ]y 'dx'] 


The second integral of the preceding expression can be evaluated once and for all at the begin- 
ning of calculations since it depends upon the airfoil boundary conditions. For the oscillating flap, 
the unsteady boundary conditions are the same above and below the airfoil, and hence the second 
integral vanishes. Equation (112) then takes the form 


^lx + ic *^l 


jdx + Atpj.i//y 


« 0/4^5 (xy f l 
+f [(7 + l)X x <P0x'^lx' +icJ (7- ^oxx'^!*^^} 


(114) 


Because the wake integral in equation (110) is slowly convergent, it is best to formulate the 
mesh boundary conditions in terms of the function + ico<p| . Let i = 1 be the column of x for 
the upstream boundary. Since derivative boundary conditions are best formulated between mesh 
points we have at x = (xj + X 2)/2 and all yj along the upstream boundary the following difference 
formula for the pressure function: 


/ X 1 + x 2 \ 

•Plx + iw*! = Op 2 j - *>ij)/( x 2 ‘ x l> + M*2j + ^lj)/ 2 = p ( — f y j/ P lj ( 115 > 

where P(x,y) represents the right-hand side of equation (114). 
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Solving for <pj • yields 

^lj = c kl^2j' c k2 P lj 

where 

c ki = (1 + iw5i/2)/(l - ico5 j/2) 
c k2 “ ^ i /C 1 ■ j/2) 

and 

*1 =x 2‘ x l- 

Substituting ip.j • into equation (84) for i = 2 yields 

a j*2j-l - (aj + bj + Ej + E 2 - E 2 c kl - q ij /2V 2j + b^ 2j+1 

= -Ei^3j + E 2 c k2 p ij 

Similarly for i = i max , the downstream mesh boundary, we have 


(116) 


(117) 


(118) 


^maxJ ^max"^ 

x] r xj ~ 

max max 


P'\ i + 'Pi .ii / X; + x 

. hnax J 'max 1J _ „( 'max 

10J 2 P V T 


max 


^ *maxJ 


(119) 


Solving for & ; yields 

'max J 

^maxj C ^^^max‘^ ^4 *maxJ 

where 

c k3 = (l -iw6 2 /2)/(l + ico5 2 /2) 

c^4 = §2/(1 + 1^082/2) 

* 

5 2 = x i ‘ x i -1- 
* 'max 'max 1 


( 120 ) 


70 



Equation (84) for i = i m ,, v - 1 with y-. : replaced by equation (120) leads to 

IIIaA ‘max- 1 

-11-1 -( a i + b j + Ei +E 2 -Ei c k3-qi -li/2Vj _ij 
J ‘max J 1 J J z 1 ‘max 1J max 1J 

(121) 

+ Vi max -ljH = ' E > C M P imaxi ' E ZW» ' 

Equations (118) and (121) follow the same pattern as equation (84) for interior columns if we 
define in equation (84) for the right-hand sides 


*lj c k2 P lj 


'Pi i ^k4 P i i 
‘max- 1 ‘max J 

and add the terms ^2^kl an< ^ p 3^"k3’ res P e ctively, to the coefficients of the diagonal terms. 

The boundary conditions on the upper and lower mesh boundaries must be given in terms of 
</? j . If the value of *p j at one point is given on the lower boundary, then the values of tpj at other 
points on the lower boundary may be found by integrating 

<Pi x + icoi £i = p ( x > y) 


where P(x,y) is defined by the right-hand side of equation (114). Between xj.j and x^, integration 
yields 



e icox p( x ' j yj)dx'. 


When the integral is evaluated by the trapezoidal rule at the Xj mesh points, then a simple recursion 
formula results from ipjj in terms of ^j-n and P(x,y). Thus 

-iOO(X: - X: i ) -ito(X; - X: i ) 

<Pil =( Pi-H e + [P(x i ,y 1 ) + P(x i . 1 y 1 )e 1-1 ] (xj - x i . 1 )/2. (122) 

For backward integration, we can solve for once ipjj is known by 


v’mi - m e 


ioj(Xj - Xj.]) 


- [PCxj.yj) e 


ico(Xj - Xj.j) 


+ P(x i . 1 ,y 1 )](x i -x i . 1 )/2. 


(123) 


Similar equations may be written for the upper boundary at y = y: 

J max 
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Since the value of must be given at one point on each of the upper and lower boundaries, 
it is wise to choose the points at which the wake integration can most easily be performed. The 
wake integral is 


(iA<p t /4s/K)(9/dy)/* e + X ^ Hq( 2 )(XjV(x - x') 2 + Ky 2 )dx'. 

J \ 


We note that for x = 1 , the integral can be simplified to 

(9/9y /~ g-iwz/O - M 2 ) h 0 ( 2 )(X 1 V z 2 + Ky 2 )dz. 

•'o 

This integral converges like 1/z 2 , but it can be expressed as an integral with exponential convergence 
by using contour integration. Consider the contour in the complex plane of figure 20 for R going 
to infinity. Since no poles lie inside the contour of C = + Cj + C 2 , then 



Figure 20. -Contours of integration in the complex z plane 
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Since the integral over Cr vanishes as R goes to infinity, then we have 


f e -icoz/(l-M ) "^/z 2 + Ky 2 )dz = -i f ^ ' ^ 2 ) Hq( 2 )(\ j ^Ky 2 - t 2 )dt 

•'0 * / 0 

+ (2/ir)/ e' w VO ' m2 ) K 0 (X, Vt 2 - Ky 2 )dt. 

•//Fi, 


</Ky 

Differentiation with respect to y yields 


(3/dy )/ e* iwz /(l 

* / 0 


o ( %,^ 


z 2 + Ky 2 )dz 


= lim { - (9 /9y) 


6 0 


i . m2) h (2) (Xl ^Ky2, t 2)dt 

• / o 


+ (2/tt)(9 /9y)/* e' wt /( 1 - M 2 ) k 0 (Xj V t 2 - Ky 2 )dt } 

JyKy+e 

\,n/k yjif 7r/2 e-^VK ycos0/0 2 Hi (2) (X jK ysin0)d0 

*^0 

+ (2/7 r)/ e' CJ '^y cosh0 ^ 2 K 1 (Xv/K’ ysinh<9)d0| 

*'0 f 


where 


j3 2 = 1 -M 2 . 


It is convenient to combine the real and imaginary parts in the form 

(9/9y )/ e‘ iojz / ^ H Q (2 )(X j v/z 2 + Ky 2 )dz 
•*0 

Xj nAk y j i ft e‘ w Vk" y cos 0/ /3 2 jjfXj y s i n d)dd 
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+ / 2 [e' w ^ ycosfl/0 2 Yi (Jk y Xj sin 0) + * cosh ^ KjCXj Jk y sinh 0)] d0 

■'0 

+[ y cosh ^ 2 K,(Xs/K ysinh0)d0| 

4/2 1 

All the integrals are defined on the range of integration, except the second integral at 9 = 0. The 
quantity in square brackets, however, goes to zero like d log 6 as 6 goes to zero. Assigning the 
value of 0 for the integrand at 6 = 0, then defines the integrands for all 0 in the range of integration. 

Analysis of Truncation Terms in the Difference Equations 

Whenever a differential equation is replaced by a difference equation and some iteration scheme 
is employed, a question arises concerning the stability of the calculations and correctness of the final 
solution. When the quantities in the difference equations are expanded in Taylor’s series about the 
central point and all terms proportional to powers of the mesh spacing are neglected, the original 
differential equation is obtained. A higher order differential equation results, however, when the 
lowest order terms in the mesh spacing are retained, and the nature of this equation will yield some 
information on the stability of the numerical method. To investigate the method set forth in this 
document, we consider the differential equation for the harmonically oscillating airfoil and will 
confine our study to an equally spaced x grid apd an equally spaced y grid, but Ax not neces- 
sarily equal to Ay. For the elliptic points, we obtain for the point i j 

u i+‘/ 2 j ( ^i+lj - *>ij)/2Ax 2 - u p/2 j (v>y - *>j_ij)/2Ax 2 

- (ioj/e)[<p i+1 j /2Ax (124) 

+ ] /2Ay 2 + (q i j/2)<^j = 0. 

We now expand in Taylor series about the point i j, letting ujj = u and <pj: = <£j . Substituting 

*i+lj = *1 + Ax ^lx + (Ax 2 /2Vi xx + ' ‘ ’ 

*>i-lj = *1 ' Ax ^lx + (Ax 2 / 2 )i£ixx + * * * 
u i+ i/ 2 j = u + (Ax/2)U X + (Ax 2 /8)u xx + • • • 
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u i-i/ 2 j = u - (Ax/2)u x + (A x 2 /8)u xx + * * * 

V?ij-1 = V\ - Ay^ ly + (Ay 2 /2Viyy + * * ' 

V’ij+l = + Ay^iy + (Ay 2 /2)v? lyy + * * • 

into equation ( 1 24) and retaining lower order terms in Ax and Ay yields 
(3/9x)(uy? lx ) - 2ico^j Je + yy + q^j 

+ (Ax 2 /12)[u^ lxxxx -2u x ^ lxxx -4i^ lxxx /e + 3u xxl p lxx /2 + u xxx *> lx /2] 

+ (Ay 2 /12) l p lyyyy = 0 

To determine whether the differential equation is elliptic or hyperbolic depends only upon the 
highest order differential operator; this is seen to be 

Ax ll( £lxxxx + ^y V 1 yyyy , 

Since u is greater than 0, the differential equation is elliptic and hence there are no characteristic 
directions to place restrictions on Ax and Ay for stability. * 

For hyperbolic points the difference equations become 

u i-l/2j^ij - y>i-ij)/2Ax 2 - Ui_3/2j(*j +1 j -^i. 2 j)/2Ax 2 

- (ico/e)[3( v ? ij - ^.jj)/2Ax - <v Mj - ^. 2j )/2Ax] 

+ [^j.j - 2v?ij + <p ij+1 ] /2 Ay 2 + qy^j/2 = 0. 

Expanding about the point i j yields the differential equation 

(uv?lx) x - 2i ^ix/ e + ^lyy + q ^l 

- Ax {ui£j xxx + lower order derivatives in x} = 0 
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The lower order terms omitted have coefficients which are negligible in comparison with the corre- 
sponding coefficients of the same derivative terms in the differential, equation. Sichel (ref. 17) 
studied the viscous transonic equation 

^x^xx ' ^yy " v ^\x\ 

and investigated the nature of shock waves. Since u is less than 0, the dominant terms of the dif- 
ferential equation become 

(-uV] xx " yy = (' u )Ax<£> j xxx + lower order derivatives 

Thus, the term (-u)Axi£j xxx introduces a numerical viscosity which Murman (ref. 2) found caused 
shock waves to develop smoothly in the iteration. From the foregoing analysis, the chosen method 
of differencing appears to be a stable formulation. 


Solution of the Complete Set of Difference Equations 

It was pointed out earlier in the text that the complete set of difference equations are linear in 
the <^jj and tfiat consideration of solving this complete set instead of column relaxation might be 
feasible. The resulting matrix is exceedingly large since there is an equation corresponding to every 
interior point of the rectangular mesh. For the computed cases, this is 73 x 56 = 4088 variables and 
equations. However, special methods for solving this set may be sought which take advantage of 
the particular sparsity of the matrix. To study the form of the equations we see that, for the 
interior points of the mesh, the difference equations for elliptic points may be written 

a j^ij-l + c ij^ij + Vij+1 + D ij^i+lj + E ij^i-lj = R ij (125) 

and for hyperbolic points 

a j^ij-l + c ij^ij + Vij+1 + E ij^-lj + F ij^i-2j = R ij- (126) 

The indices i j here run from 1 to i max and j max for interior points of the mesh. 

For 0 values of the subscript, the terms vanish on the left side and become incorporated into 
the Rjj since they are mesh far-field boundary conditions. Since there are i max columns and j max 
rows in the interior of the mesh, then there are N = i max x j max variables with N equations to be 
solved for them. 
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We order the variables as follows: Let 


X n ~'P\) 


where n = i + 0 - l)i max with 1 <i<i maY , 1 <j <j 


max’ 


J max - 


Then equation (125) takes the form 

a j X n-i max + E ij X n-l + c ij X n + D ij X n+l + b j X n+i max = R ij 
for the n* b equation at elliptic points, and equation (126) takes the form 


(127) 


a j X n-i max + F ij X n-2 + E ij X n-l + C ij X n + b j X n+i max R ij 


(128) 


for the n* b equation at hyperbolic points. Since only the interior points about the airfoil are apt to 
be hyperbolic, we will study the form of the equations assuming all elliptic points. Equation ( 1 27) 
will contain all five terms on the left side when 


•max < N - i max 

For the other equations, we will have one negative subscript for X which represents the mesh far-field 
boundary conditions and will be incorporated into the right-hand side term, Ry. All main diagonal 
terms are nonzero, and we see that the coefficients of the terms adjacent to the diagonal terms are 
also nonzero. Also, the term i max points to the left of the diagonal and the term i max points to the 
right have nonzero coefficients when 


•max < n < N - i max 
Thus, the matrix has only five diagonals. 

As seen from equation (128), the addition of a supersonic point in the flow field places a zero in 
the diagonal to the right of the main diagonal and adds a nonzero value in the second diagonal to the 
left of the main diagonal. Thus, a local supersonic region introduces a sixth nonzero diagonal to the 
matrix. 
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